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ABSTRACT 

We investigate shear and buoyancy instabilities in radially stratified, magnetized, 
cylindrical flows, for application to magnetocentrifugally driven winds - such as those 
from protostars - and to magnetized accretion disks. Our motivation is to characterize 
the susceptibility of cold MHD disk winds to growing internal perturbations, and to 
understand the relation of wind instabilities to known accretion disk instabilities. Using 
four different linear analysis techniques, we identify and study nine principal types 
of unstable or overstable disturbances, providing numerical and analytic solutions for 
growth rates for a wide range of parameters. 

When magnetic fields are predominantly toroidal, as in protostellar winds far from 
their source, we find the system is susceptible to growth of five different kinds of pertur- 
bations: axisymmetric fundamental and toroidal resonance modes, axisymmetric and 
non-axisymmetric toroidal buoyancy modes, and non-axisymmetric magnetorotational 
modes. Winds having a sufficiently steep field gradient (din B/ din R < —0.75 for a 
purely toroidal-field case) are globally unstable to the long wavelength fundamental 
mode concentrated at small radii; these promote the establishment of narrow dense jets 
in the centers of wider winds. Long-wavelength outer-wind modes are all stable for 
power-law wind equilibria. The toroidal buoyancy instabilities promote small-scale ra- 
dial mixing provided the equilibrium has nonzero magnetic forces. For low-temperature 
toroidal-B winds, both axisymmetric and non-axisymmetric magnetorotational insta- 
bilities have very low growth rates. The stabilization of buoyancy instabilities by shear 
and of magnetorotational instabilities by compressibility may be important in allowing 
cold MHD winds to propagate over vast distances in space. 

When magnetic fields are predominantly poloidal, as may occur in protostellar winds 
close to their source or in astrophysical disks, we find the system is susceptible to 
four additional growing modes: axisymmetric magnetorotational (Balbus-Hawley) , ax- 
isymmetric poloidal buoyancy, non-axisymmetric geometric buoyancy, and poloidal res- 
onance modes. The well-known axisymmetric Balbus-Hawley mode has the fastest 
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growth rate. When the magnetic field is nonuniform, the axisymmetric poloidal buoy- 
ancy mode promotes radial mixing on small scales. The geometric poloidal buoyancy 
mode requires high m, thus is readily stabilized by shear. 

Previous work on magnetorotational instabilities has concentrated on near- 
incompressible systems (accretion disks or stellar interiors). We extend this analysis 
to allow for compressibility (important in winds). We introduce a "coherent wavelet" 
technique (a WKB temporal approximation), and derive closed- form analytic expres- 
sions for instantaneous instability criteria, growth rates, and net amplification factors for 
generalized non-axisymmctric magnetorotational instabilities in compressible flows with 
both poloidal and toroidal fields. We confirm that these are in excellent agreement with 
the results of shearing-sheet temporal integrations, and that "locally-axisymmetric" 
perturbations have the largest amplifications only provided (k • va)/J^ ^ 1- 

Subject headings: accretion, accretion disks — ISM: dynamics — ISM: jets and outflows 
— ISM: magnetic fields — MHD — stars: pre-main-sequence 

1. Introduction 

The ubiquity of energetic molecular outflows and atomic jets from young stellar objects (YSOs) 
ranging from deeply embedded infrared sources to classical T Tauri stars suggests that they are an 
inescapable by-product of star formation (e.g.. Richer et al. 2000, and references therein). Winds 
from YSO disks play an important role in shedding angular momentum carried by inflowing mate- 
rial, thereby permitting further accretion in order for the central objects to attain stellar dimensions 
(Hartmann & MacGregor 1982; Pudritz &: Norman 1986; Shu et al. 1988). These winds may also 
have strong effects on the dynamical evolution of the parent cloud by providing a source of turbu- 
lent energy (Norman & Silk 1980), and may help to determine the final masses of stars by reversing 
the infall of surrounding gas (Shu, Adams, & Lizano 1987). Therefore, understanding the physics 
of protostellar winds is of essential importance to the theory of star formation. 

Among the various scenarios regarding the origin and nature of the protostellar winds, the 
most promising is magnetohydrodynamic (MHD) models in which winds are driven by the inter- 
action of the centrifugal force with open magnetic fields threading rapidly rotating disks. These 
magnetocentrifugally driven wind models can account for observed high mass and momentum losses 
(Lada 1985; Bachiller 1996) and the kinematic and structural characteristics of bipolar molecular 
outflows in general (Li &: Shu 1996), as well as in specific cases (e.g., HHlll, cf. Nagar et al. 1997). 
It is, however, still controversial whether the wind originates only from a small magnetosphere-disk 
interaction region near the central star (Shu et al. 1994, 2000), or whether it emanates from an 
extended region of the disk (following the seminal model of Blandford & Payne 1982; see, e.g., 
Konigl & Pudritz 2000). 

Although the role of magnetic fields in driving protostellar winds is by now well established 
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(at least theoretically), their complementary role in governing the stability properties of winds is 
less well explored. In addition, azimuthal and vertical shear within winds may also affect their 
stability properties. Questions of wind stability are potentially important for both large scale and 
small-scale phenomena. These include understanding the role of magnetic fields and velocity shear 
in (a) helping winds to propagate over enormous distances (up to a factor ~ 10^ in dynamic range) 
through the ISM in parsec-scale giant HH flows (Reipurth, Bally, k, Devine 1997); (b) creating bright 
HH "knots" spaced throughout optical jets (Hartigan et al. 2000); (c) governing the angular extent 
of the emergent wind and establishing the momentum distribution for driving molecular outflows; 
and (d) converting large-scale ordered flow energy to jet heating through small-scale instabilities. 
In addition, it is important to study the dynamical properties of large classes of theoretical wind 
solutions to test whether they are stable equilibria which can represent real astronomical systems, 
or whether they are unstable equilibria which should rarely be observed in nature because they 
evolve rapidly into other configurations. 

Previous studies of time-dependent behavior in steadily-input MHD winds have focussed pri- 
marily on the Kelvin-Helmholtz instabilities driven by the interaction of the wind surface layer 
with the ambient medium (see, e.g., the studies of Appl & Camenzind 1992; Hardee ct al. 1992; 
Rosen et al. 1999, and references therein). Generally, heavy jets containing strong toroidal fields 
are relatively resistant to these instabilities. In addition to "driven" instabilities resulting from 
boundary conditions, winds may also be subject to "free" instabilities in their interiors. Whether 
a given wind solution is internally unstable must depend on the velocity shear, magnetic geometry, 
and internal stratification. One route to studying internal wind instabilities is via time-dependent 
numerical simulations of winds. Although such studies (e.g., Ouycd & Pudritz 1997, and references 
therein) have yielded intriguing results on the development of episodic knots in MHD winds, the 
computational demands in carrying out simulations precludes extensive exploration of parameter 
space, large spatial dynamic range, or very long-term integration. In addition, some of the time- 
dependent internal features found in simulations may be introduced by particular choices of inflow 
boundary conditions that are inconsistent with a steady-state flow, rather than occurring as a result 
of intrinsic instability of the wind. 

Due to the importance of gaseous accretion disks in a wide variety of astrophysical systems, 
major attention has focussed on disk dynamics, and, in particular, the role of instability-driven 
turbulence in angular momentum transport (see, e.g., Balbus & Hawlcy 1998; Stone et al. 2000, 
for reviews). Saturated magnetorotational instabilities (hereafter MRIs; Balbus & Hawley 1991, 
1992; Hawley & Balbus 1991; Hawley, Gammie, & Balbus 1995) represent perhaps the most impor- 
tant local dynamical process affecting disk evolution. Magnetized disk winds share many generic 
properties with disks, so it is interesting to investigate the potential importance of MRIs in winds. 

In this work, we investigate the internal stability of rotating, magnetized protostellar winds 
to (primarily) local shear and buoyancy modes. We also extend previous studies of local MRIs in 
accretion disks. The fundamental difference between "wind" and "disk" systems in our idealized 
models is in the absence or presence of gravity as a conflning force. These systems may also be 
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distinguished by the geometry of the magnetic field, with toroidally-dominant fields expected in 

the wind case, but cither poloidal or toroidal fields possible for the disk case. Our most general 
analysis and results apply to cold flows, but we also perform separate calculations (see §7) including 
thermal effects which specialize to local analysis of MRIs. 

Whether a wind originates from a narrow boundary layer or an extended radial region, the 
radial expansion of the flow will lead to shear in the azimuthal velocity fleld of the asymptotic state. 
The total specific angular momentum of the flow, J = R{vip — BpoiB^/{4TrpVpoi)), is conserved along 
streamlines (where -Bpoi, and Vpoi, v^f, are the poloidal, toroidal components of the magnetic fleld 
and the flow velocity). If the Alfven mach number Ma = Wpoi/^A,poi S> 1, the kinetic part dominates 
the speciflc angular momentum and v^f) ~ J/R. Thus, the angular velocity Q = v^/ R J/R? in 
the asymptotic wind will have a gradient dlnVl/ d\\iR = din J /din R — 2. If the wind comes 
from a boundary layer, then din J / din R may be small; if the wind originates from a large region 
with self-similar scalings, then d\nJ/d\nR = 1/2. In either case, dlnQ,/dlnR is expected to be 
a negative, order- unity quantity for wind systems. For thin disk systems (i.e., negligible pressure 
support), a Keplerian radial proflle d\nfl/d\nR = —3/2 is expected if the central mass is dominant. 
For the analysis of this paper (except where noted otherwise), we adopt din 0,/ din R = —3/2 for 
both "wind" and "disk" systems, but our qualitative results are insensitive to this assumption. As 
discussed below, order-unity radial logarithmic gradients may also be expected in the magnetic 
fleld strengths; we allow for a range of magnetic gradients. 

Signiflcant shear may also exist in the poloidal velocities of jets if they originate from an ex- 
tended radial region. The asymptotic outflow speed vqz generally scales linearly with the rotational 
speed at the footpoint -Rfoot of the streamline, so that 

dlnvpz _ dlavQ^ ^ln^;foot Olnfifopt ^ 1 dlnRfpot 
dlnR din Ufoot 9 In i?foot dlnR 2 dlnR 

If the range of footpoint radii is small compared to the range of asymptotic radii (as for a wind from a 
boundary layer), then \dlnvQz/dlnR\ <C 1; if the radial ranges arc comparable, then 51nfoz/91ni? 
is negative and order unity. For disks, the vertical velocity shear is negligible. We allow for a range 
of vertical shear rates in the present analysis. 

Our analysis consists of developing and solving sets of linearized MHD equations for several 
general classes of background flows. Both axisymmetric and non-axisymmetric disturbances are 
explored. Even in linearized form, MHD problems present considerable technical challenges. Thus, 
instead of attacking sophisticated sets of steady-state wind solutions, for many speciflc examples 
we will take as an unperturbed configuration one of the power-law cylindrical equilibrium solutions 
recently identified by Ostriker (1997) as asymptotic states of self-similar disk winds. These have 
density p oc R~'^, B oc v oc and sound speed Cg = 0. Although these adopted 

initial configurations are relatively simple, they retain general asymptotic characteristics of MHD 
disk winds in the sense that they have both azimuthal and vertical magnetic field and velocity 
components with arbitrary ratios, and all the physical variables have radial gradients. These 
gradients may also be thought of as representing local scalings within a more complex overall 
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stratification. We also include models without vertical motion but with significant equilibrium 

gravity to study stability in magnetized astrophysical disks. In our analysis of MRIs, we use 
equilibria with uniform B,p, and Cg 7^ 0, and take O, oc with arbitrary a. 

Because the systems we are studying contain significant azimuthal shear, an arbitrary initial 
spatial planform is not maintained indefinitely. When Q' ^ and the azimuthal wavcnumber m ^ 
(and/or when v'q^ 7^ and the vertical wavenumber kz 7^ 0; the prime represents a differentiation 
with respect to R) , spatial wavefunctions may become increasingly radially corrugated in time due 
to the kinematic shearing of the wavefronts imposed in the initial conditions. If we describe the 
radial spatial wavefunction in terms of the amplitudes of Fourier coefficients with radial wavenum- 
bers /cr, this corresponds at late times to a secular increase in the amplitudes of large-A;R terms and 
decrease in the amplitudes of small-A;R terms. As we shall show, all the disturbances we identify 
are stabilized at sufficiently large k^. Thus, if mO' and/or kzVQ^ 7^ 0, the net amplification factor 
for any arbitrary initial disturbance is limited by the rate of kinematic growth of radial corrugation 
compared to the growth rate of any dynamically-driven instabilities. 

In previous work, two complementary analytical methods have been used to study small- 
amplitude disturbances in shearing astrophysical systems. One approach adopts the "shearing- 
sheet" formalism, and integrates the local, time-dependent, linearized wave equations directly to 

obtain the evolutionary behavior of shearing wavelets treated as an initial- value problem (Goldreich 
& Lynden-Bell 1965; Julian & Toomre 1966; Balbus & Hawley 1992). An alternative analytical 
approach uses WKB techniques to derive dispersion relations for spatial Fourier modes, superposi- 
tions of which represent local wavefunctions (e.g., Shu 1974, 1992; Ryu &; Goodman 1992; Terquem 
&L Papaloizou 1996). This paper includes analyses using both approaches, and also introduces a 
hybrid technique which we term a "coherent wavelet" analysis. We adopt the "modal" strategy 
in order to identify characteristic instantaneous growth rates and physical mechanisms for a wide 
variety of spatial disturbances. By considering the time over which a spatial pattern is altered 
by shear, we can estimate the net amplification factor of a given initial modal disturbance. We 
use the shearing-sheet formalism for studying magnetorotational instabilities, which are cut off at 
relatively small values of Rk^^/m (whereas the modal analysis applies to large Rk^/m), and also for 
studying buoyancy instabilities in the high-m regime where modes are short-lived. We show that 
the results obtained from the shearing-sheet integrations in both cases are in excellent agreement 
with the predictions of a coherent wavelet analysis, in which time-dependent growth rates ^(t) 
are computed by time-localizing the shearing-sheet equations and solving an analytic dispersion 
relation. 

The organization of this paper is as follows: We begin by studying instabilities in cold, magne- 
tized winds. In §2, the basic MHD equations and the specific adopted wind equilibria are described. 
In §3, we analyze the stability of winds to the simplest perturbation with fcz = m = 0, where kz and 
m are respectively the vertical and azimuthal wavenunibers of the perturbation. We term these 
the "fundamental modes" ; we present solutions for stable and unstable global modes under the as- 
sumption of free Lagrangian boundary conditions. The modal analysis and general local dispersion 
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relation for cold flows with arbitrary and m are presented in §4. We present numerical solutions 

of the dispersion relation for both axisymmetric and non-axisymmctric perturbations in §5. In §6, 
we classify the unstable or overstable modes and provide the physical interpretation for each mode. 
Next, we include (variable) thermal pressure terms to compare the susceptibility of cold winds vs. 
warm disks to shear instabilities. In §7 we analyze the axisymmetric Balbus-Hawley instability of 
poloidal fields and the non-axisymmetric MRI of toroidal fields, discuss the respective instability 
mechanisms, and provide the corresponding instability criteria. Here, wc use the coherent wavelet 
technique to compute growth rates, and compare with direct shearing-sheet integrations. The gen- 
eralized instability criteria and net amplification factors for the magnetorotational disturbances 
with both toroidal and poloidal background fields are also derived. Finally in §8, we summarize 
and discuss conclusions of the present work. 



2. Basic Equations and Cylindrical Equilibrium for Cold Wind 

We begin with the ideal MHD equations 

§^ + V-(pv)=0, (1) 

dv 1 VP 

+ vVv = — (VxB) xB V$G, (2) 

at Airp p 

— = Vx(vxB), (3) 



and 



V • B = 0, (4) 



where p is the density, v is the fluid velocity, B is the magnetic field, P is the thermal pressure, 
and — V$G = — g is the gravitational force due to a central object. We ignore self-gravity in the 
flow. 

We now consider cold, magnetized cylindrical flows. Since the flow velocity in disk winds is 
always supersonic except in the vicinity of the disk where the material is lifted by the thermal 
pressure (Blandford & Payne 1982), the thermal pressure term in eq. (2) can generally be neglected 
compared to magnetic stress. Except for investigations of generalized MRIs (§7), we shall drop the 
thermal pressure term. We adopt standard cylindrical coordinates {R, (p, z) . 

By assuming that = B^i = and all variables are independent of z, we have a general 
equilibrium condition from eq. (2) 

vl 1 (Bl \ 

where a prime denotes differentiation with respect to R. At a large distance from the origin, the 
gravitational force due to the central source can also be ignored on the grounds that magnetic and 
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centrifugal forces far exceed it. In this case, the magnetic hoop stress acting inward is the only force 
that balances the outward centrifugal force and outward magnetic pressure gradient force (under 
the assumption that the magnetic field strength decreases outward). 

As an initial equilibrium configuration of the wind, for specific cases we will adopt the asymp- 
totic solutions for cylindrically symmetric axial flows presented by Ostriker (1997). All variables 
have power-law dependences on R: p oc R~'i,B^ on on i?"^^"'"'?)/^, and oc fz oc R~^l^. We 
define the local pitch angle % of the magnetic fields as i = tan~^ {B^jB^. Neglecting in eq. (5), 
radial momentum balance requires 

^1 = Y(cos2i-g), (6) 
where is the local Alfven speed defined by 

I'A = + ^Az with = ■^^= and vaz = 

It is obvious from eq. (6) that there would be no such power-law solutions q > 1. This is 
because when g > 1, the magnetic field has so steep a gradient that the corresponding pressure 
force always exceeds the tension. Therefore, to ensure force balance and cylindrical collimation in 
winds with power-law profiles, the magnetic field strength must decline with R more slowly than 
R-\ 

We can define an angular velocity Jlf = ^2 — VpoiB^/{RBpQi) as that of a rotating frame in 
which the fiow of winds is parallel to the local field line, ilf is the rotation rate of the magnetic 
field lines thought of as rigid wires. In such a frame, the family of solutions can be completely 
described in terms of scaled values of the specific angular momentum j, the Bernoulli constant e, 
and q, where 

(Ostriker 1997). The condition for a super- Alfvenic outflow velocity requires < j < 1. Generally 
speaking, the pitch angle i does not depend on q, although one can parameterize i in terms of q, e, 
and j. However, for flows originating from a Kepler-rotating disk, angular momentum and energy 

conservation requirements limit the range of i available to an equilibrium (asymptotic) magnetic 
field configuration. Utilizing eqs. (15) to (23) of Ostriker (1997) one can show that the maximum, 
over all permitted values of e and j, pitch angle imax is given by 

tan^ imax = f j -1, (7) 



which is attained when e = and j = 1- If j > 1, the streamline never reaches the Alfven radius. 
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3. Fundamental Mode 

3.1. Dynamical Equations 

We first consider the response of tfie equilibrium state when small, axisymmetric perturbations 
with an infinite wavelength along a vertical direction are imposed. We term the waves with = 
and m = the "fundamental modes" analogous to eigenfunctions of oscillations without any node. 
Let the subscripts and 1 denote the equilibrium and perturbed states, respectively. Linearizing 
the set of the dynamical equations (1) to (4), we may write 



dv 



IR 



dt 



= 2nvi^ - 



47rpo 



dpi ^9 



dv 



dt 

at 



OB 



dt 



dBi, 

at 



_d_ 
~'dR 



(So^f ir) , 
{RBq^vik) , 



and Sir = 0. In eq. (10), k stands for the epicycle frequency 



R2 




+ Bo • Bq 



(8) 
(9) 

(10) 

(11) 
(12) 

(13) 



Combining the perturbed equations (8) — (13) and eliminating all other variables in favor of 
the perturbed radial velocity f ir, one obtains the wave equation 

1 a^viK 



vl at^ 



a-'vi^ ^dln{RB^)dviK 



aR^ 



dR 



aR 



1 



i?2 



-Q- + 1 + cos i 



For power-law profiles, this becomes 

1 a'^viK a'^viK q dviR 



af^ 



R dR 



.dlnpo 

d\nR ) RBl dR \pQ 



« 1/1-9 2. 



ViK- 



(14a) 



(14b) 



To transform eqs. (14) to the Schrodinger form, we define a new independent variable ^ through 



or vvR = R'i''^-^{R)e 



iujt 



(15) 
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for the power-law case. Then, we have 



with K{R) defined by 



+ ii:2(^)^ = 0, (16) 



or 

for the power-law case. Local (WKB) solutions to eq. (16) have * ~ e^fkndR ^j^j^ j^f^^ i g^^^j 
dk^/dR ^ A;^. In this case, ^" —k^^, and we can use K^{R) = k^ to write a local dispersion 
relation 

= vikl, (18) 

which corresponds to MHD fast modes propagating along the radial direction. When is com- 
parable to or smaller than v\/ E? , however, modes are not localized, and solutions must be sought 
as a global problem subject to boundary conditions. 



3.2. Global Analysis for the Fundamental Modes 

In the previous section we showed that there is no short wavelength (local) unstable funda- 
mental mode with fc^ = m = in self-similar MHD disk winds. Here, we present the results of a 
global normal-mode analysis performed with carefully chosen boundary conditions, and adopting 
the power-law equilibrium. Define the dimensionlcss radial variable r = R/Re, and dimensionless 
parameters a = (q — 1) cos^ i + q{2 — q)/A and = uj"^ R^ / v\{Rc) , with R^ being the position of 
the unperturbed outer edge of the wind. Then eq. (16) can be cast into the form 

72, 



d^'l! ( a 



, a2r)* = 0. (19) 
It is not difficult to show that the general solutions of eq. (19) are 



^ = 1 



' A^J^{2ar^''^/i) + B^^{2ar^/^/i), if > 0, (20a) 
C^h{2\u\r^l^/^) + D^K^{2\u\r^/'^/Z), if < q, (20b) 



with "iu = Vl ~ 4q; = y^(l — qY + 4(1 — q) cos^ i. In eqs. (20), J^-, Yy, and 1^, are the ordinary 
and modified Bessel functions of the 1st and 2nd kinds, respectively, and the coefficients A, B, C, 
and D are constants to be determined from imposed boundary conditions. 
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Let us consider the case of a free Lagrangian boundary at which the total pressure due to 
initial and perturbed fields balances with a fixed external pressure at both inner and outer edges, 
which is equal to the unperturbed magnetic pressure. If the total pressure at an edge of an outflow 
is different from the external pressure, the boundary itself will move until a new balance exists. To 
first order, this condition of constant pressure at the boundary is written 

where i?b is the location of the perturbed boundary. The first term of eq. (21) represents the change 
in the total pressure due to the boundary movement, while the second term arises from perturbed 
magnetic pressure itself. All quantities are evaluated at the unperturbed boundary position. Using 
eqs. (12), (13), and (15), and using dR\)/dt = uir at boundaries, we find the desired boundary 
conditions are 

all + sin^ i ^ , , , 

— + ^ = 0, at r = n and 1, (22) 

dr r 

where r\ = Ri/R^ is the normalized distance of an inner boundary from the axis. Together with the 
boundary conditions (22), eq. (19) forms a Sturm-Liouville system. By employing the variational 
principle one can show that cr^ is real and that <t^(^') is stationary subject to an arbitrary variation 
of ^. 

When (T^ > (stable modes), the oscillatory properties of Jy and Yy guarantee the existence 
of discrete eigenvalues with n denoting the number of nodes in the corresponding eigenfunction 
*n- The resulting eigenvalues for n = 10~^ and 10~^, and < z < imax(<?) are plotted in Fig. 1. 
Only a few cases with small n are shown. Eigenvalues associated with different Vs fill each shaded 
area completely. When ri = 10~^, eigenfunctions which link inner and outer boundaries have to 
extend across enormous changes in density and magnetic field strengths. In this case, B/A -^1 
and eigenvalues become rather insensitive to the local properties such as q and i. When r\ = 10"^, 
however, the wind mimics a slender hollow cylinder. The variation in density and field strengths 
over radius is slight, causing eigenvalues to be sensitive to i and q. In addition, the narrow width 
of the wind changes the number of nodes in eigenfunctions. When q > 0.5, for example, the 
eigenfunctions with ri = 10~^ have almost the same eigenvalues as, but one more node than, the 
n = 10~^ case, as seen in Fig. 1. 

When n ^ 1 and cr^ 3> 1, the asymptotic sohitions to eqs. (20a) and (22) gives = 37rn/2 + 
2>tt{2u + l)/8. These are plotted with dotted lines in Fig. lb, and show good agreement with the 
values calculated without any assumption (even for n = 0). The case with q = 0.5 and i = is 
special, because the slope of the eigenfunction at the inner boundary is —1/4, which automatically 
satisfies the boundary condition (22). In this case the asymptotic eigenvalues are = 37rn/2, 
drawn as filled circles in Fig. lb. Eigenvalues have no upper bound as n — oo, which is a general 
property of solutions to a Sturm-Liouville equation (Morse &; Feshbach 1953). 

Now consider the unstable global solutions with cr^ < 0. Let ■^i and V'2 be the two linearly 
independent solutions of *: Vi = ^A^i^l/(2|<J|r^/2/3) and V2 = v^-fCi/(2|cr|r^/^/3) such that * = 
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Cipi + Dip2- Because %l>i is a monotonically increasing function of r (i.e., Vi^V'i > always) and 
^ must have a negative logarithmic slope at the inner and outer boundaries (cf. cq. [22]), tpi alone 
can not constitute eigcnfunctions for global modes. In addition, since tpi increases exponentially 
for a large value of |c7|r^/^, while 1/^2; ~^ 0; outer free Lagrangian boundary condition requires 
C/D — 0. Although C is not strictly zero when |cr| has a relatively small value, the contribution 
of ipi to global solutions near the inner boundary is negligibly small. Thus, unstable eigenvalues, 
if they exist, are essentially determined by the inner boundary condition imposed on ■02- 



As we move inward from the outer boundary, tp2 rapidly increases asymptoting to 



2u 



1 _LU ^3i. , 

32'^sin(7rz/)r2(i/ + 1) 



(23) 



for r <C 1 (cf. Abramowitz &; Stegun 1965), where r(z/+l) is a Gamma function. In fact, (1— 3z^)/2 is 
the maximum logarithmic slope ip2i'r') can ever attain. Prom the inner boundary constraint (22), the 
existence of unstable global solutions is guaranteed if (1 — 3u)/2 > —{q/2 + sin^ i), or, equivalently 

,>l-^i±f^, (24) 

is satisfied. Eq. (24) is the global instability criterion for the fundamental modes of self-similar, 
cold, magnetized winds, subject to the free boundary conditions expressed by eq. (22). By putting 
C = and neglecting higher order terms in -02, we derive from eqs. (22) and (23) the approximate, 
analytic expression for the eigenvalues of the global instability 



I I 3/2 i^^l^i 

(7 T; — 



VA{Ri) 



sin{TTu)T'^{iy + 1) (1 - 3u + q + 2sin^ i) 



l/2v 

(25) 



which again shows that \a\ has a positive real value if cq. (24) holds. 



In Fig. 2 wc plot the approximate growth rates for unstable modes from eq. (25) as dotted 
lines, as well as the exact growth rates numerically computed for r\ = 10^^ (thin solid lines) and 
for Ti = 0.1 (dashed lines). The curves shown are for i = 0°,5°,- • -,35", 40" from right to left, 
and the uppermost thick lines are for Zmax calculated from eq. (7). Note that varying the width of 
outflow via Ti yields very little change in the plotted solutions: ri-dependence of the growth rates 
appears mainly through the product \(j\r-^ . Eq. (25) gives accurate growth rates for relatively 
small values of |(T|rj^''^, while its estimates deviate up to 16% from the exact values as |o"|rj^^^ 
becomes larger. In this case we need to include next order terms in ^2 (cf. eq. [23]) to obtain 
more accurate results. For a given set of equilibrium parameters, we note that whereas there 
exist an infinite set of stable eigenmodes, there is (at most) a unique unstable eigenmode. Noting 
|a;| = \a\v\{Re)/ Re = \a\rf^'^VA{Ri)/Ri, we expect from Fig. 2 that the system is typically globally 
unstable within ~5 crossing times of Alfven waves at the inner boundary. 

Once the ratios of the coefficients and the eigenfrequencies are found for the fundamental 
modes, one can easily construct radial solutions for the perturbed variables: pi, vir, vitf,, and Bi^p. 
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These are plotted in Fig. 3 for z = 0° and n = 10 ^. Fig. 3a corresponds to a stable case with 

3/2 

q = 0.4 and cjo = 2.42, while Fig. 3b depicts an unstable case with q = 0.6 and |cr|rj = 0.11. 
Although normalization is arbitrary, we note that for the unstable modes, the negative radial 
velocity case drives the entire system into a more stable configuration (with lower magnetic energy) 
when the equilibrium magnetic field is predominantly toroidal. This can be shown as follows: 
Let 5M, SEb, and S^b denote the mass, magnetic energy, and toroidal magnetic flux per unit 
height in a local flux tube. Then we have SEb = (tt /2){S^b / SM)'^ 5M pR'^ . For a given flux tube, 
6M and S^b/SM are constant in time and SM > 0. Thus, sgn d{5EB)/dt = sgn d{pR^)/dt = 
sgn [pqR^{vi'q^/ R — dvx^/ dR)\ from the equation of continuity. If vib/R dominates Ovib/OR and 
^iR < 0, then sgn d(6EB)/dt < 0; magnetic energy decreases with time, meaning that the system 
evolves into a more stable state. Thus we scale vib/vq^ = 1 at r = 1 for Fig. 3a and vib/vq^ = — 1 
at r = Ti for Fig. 3b, respectively. Note that stable eigenfunctions have their largest amplitude 
near the outer boundary, while the inner, high density region is nearly static during oscillation. 
Unstable eigenfunctions, on the other hand, are almost zero except in the region close to the 
inner boundary. The respective inner-region vs. outer-region predominance of unstable vs. stable 
eigenfunctions reflects the respective characteristic frequencies as well: the inner-wind unstable 
modes grow at large rates comparable to Alfven frequencies in the interior, whereas outer- wind 
stable modes oscillate at low frequencies comparable to the Alfven frequencies in the exterior. 

We remark that there is no globally unstable fundamental mode when one adopts rigid bound- 
aries with ^'(r) = at both r = n and r = 1, instead of the free Lagrangian boundaries, since both 
■01 and 02 arc monotonic functions of R. If ^''(r) = is imposed at both boundaries (cf. Dubrulle 
& Knobloch 1993), however, we still have unstable fundamental modes with a different instability 
criterion^ and different growth rates. 

We discuss the significance of fundamental modes to protostellar outflows in §8.3. 



4. Local Analysis for Cold Winds 



We now consider general non-axisymmetric Eulerian perturbations with small amplitudes. 
Neglecting the effects of thermal pressure and external gravity due to a central object, we linearize 
eqs. (1)~(4) 

dpi 1 _ , Id 



dt 



^^^-(i?Po^m)-^^(Po^i,) 



d_ 
dz 



{PQVlz 



(26) 



-dT = 



47rpo 

dvi<i> _ _ 
dt ~ ~ 



2Bq^Bi^ Bo^dBiB , D dBiB d x , Pi 

5 ~5 ajT + -°02;— 5 aBv^O • til H 

R R ocp dz oR po 



^°'^+Bo-B' 



R 



if_ 1 



Bib. + B[ 



Oz 



dB 



1 dBi^ 
dz R~d^ 



14, 



0- ^oj 
(27)" 
(28) 



^In this case, eq. (24) would become (1 — 3i')/4 = a > 0, corresponding to B?{K^{R) — of^ /v\) > 0. 
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dv 



Iz 



dt 



1 



dt 



B04, dviYi 



R 



+ B, 



Oz- 



1 dBu 
R d(j) 

dviR 
dz ' 



dB 



dz 



dB 



dt 
dB^ 

dt 



= Rn'BiR - ^{Bo^vm) + 5oz^ - B 



where the Lagrangian time derivative is denoted by 

d 
di 



d d d 



(29) 
(30) 
(31) 

(32) 

(33) 



and again the subscripts and 1 indicate the equiUbrium and perturbation variables, respectively. 

Since all the coefficients of the perturbed variables in eqs. (26)~(32) do not depend on the co- 
ordinates (f) and we may look for solutions having sinusoidal dependence on ^ and z. Furthermore, 
if there exist any normal modes, we can write eigenfunctions in the form 



Xi(i?,<^,^,t)=Xi(i?)e^("*'^+'^^^-"*\ 



(34) 



where xi refers to any physical variable of perturbations. Substituting eq. (34) into the set of 
eqs. (26)~(32) and eliminating all other variables in terms of the radial Lagrangian displacement 
= —viji/ioj with a Doppler shifted frequency 



00 = uj — — kzVoz, 
we obtain the second order differential equation 



where 



HiR) ^ 



with 



d^^R , d 



dR^ dR 



~,2 



dR 



~,2 



dlnpo 
dR 
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Am 

dR ^ RBl 

.2 r2 ^^A 



d 



dR 



B, 



Bo 



R 



viik'Fi:^ + ^FBG+G^ + ^ 



R 



1 



Bl 



R 



+ Bq • Bq 



' 1 ) G± = — \^Bqz ± hBo^j , 



(35) 



(36a) 



(36b) 
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oj\ = Cu'^ — v\{]s. ■ ho)"^ , = ^'^ — vpji^ , and k=(0, ^,A;z). 

R 

Here, Fb represents the equilibrium magnetic force, and and u-p are frequencies connected 
to the Alfvenic and fast magnetosonic modes in cold MHD fluids, respectively, bo = Bo/|Bo| is 
the unit vector along an equilibrium field direction, and finally k is a vector wavenumber. When 
m = = 0, only the terms in the first bracket in the definition of H{R) do not vanish, recovering 
the radial wave equation (eq. [16], [17a]) for the fundamental mode. 

The second order differential equation (35) has a singularity at u)\ = 0, but if we treat a fully 
general problem including the thermal effects of compressible gas, we will find another singularity 
(a so called cusp singularity) at the positions where Doppler shifted frequencies of traveling waves 
match with slow MHD wave frequencies of the medium (cf. Roberts 1985). For an incompressible 
medium, a;^ = singularities are often referred to as shear Alfven singularities where because of 
resonances the characteristics of waves propagating radially would be modified to be either absorbed 
into or amplified by background medium, if considered as a boundary value problem (Ross et al. 
1982; Curry &: Pudritz 1996). As pointed out by Appert et al. (1974), the locations with tDp = 
in eq. (35) are not singularities; these cut-off points in our local analyses appear as resonance 
waves with frequencies having relatively small imaginary parts, suggesting potential attenuation or 
amplification of amplitudes. 

To remove the second term in eq. (35) we further define 



^^^^![RBl%\ "\ (37) 



Then eq. (35) is reduced to the standard Schrodinger form of eq. (16), with generalized K'^{R) 
defined by 



2 



(38) 



Generally speaking, K[R) is a function of R for fixed values of m, /c^, and uj. However, we can 
still consider the behavior in a local sense near some fixed i?o, such that K is close to K{Ro). This 
is mathematically formalized as described in Lin et al. (1993) and Terquem & Papaloizou (1996). 
Let us consider in the nonuniform background the spatially localized wave packet of the form 

* = v(i? - Roy^''^^~^°'^ + 0(7^), 

where %lj{r) is a function which is non-zero only in a small neighborhood of r = i? — i?o = 0. The 
scale over which ^'(r) varies significantly must tend to zero as /cr 00, but no faster than . 
Then, to leading order, cP"^ /dR^ —k^^, and the solution k"^ = K'^{R, oj, m, k^) of the Schrodinger 
equation yields a local dispersion relation with the right hand side evaluated at a reference point Ro, 
provided /cr is limited to a sufficiently large value (i.e., i?o^R ^ 1)- We may invert this dispersion 
relation to find lo = W{m,kz,k^_,R), so that k^ now plays the role of an independent parameter 
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and the dispersion relation yields the Doppler-shifted frequency a) of a wave near a position Rq 
having local wavevector k. This is equivalent to a standard WKB approximation in the radial 
direction. 

Solution of the local dispersion relation near yields 

u = W{m, h, fcR, Ro) + m^{Ro) + hvo^{Ro) + 0{r), (39) 

where the 0{r) term is (mJ7'(i?o) + hv'Q^{Ro))r. Defining dkn/clR = -{dW /dR)/{dW /dkn) ~ 
k^i/R, the WKB condition \dkji/dR\ <C k^ will be satisfied for k^iR ^ 1. For normal mode 
solutions, the 0{r) term in uj must provide a negligible contribution to the phase; this requires 
that we must have \m^'{Ro) + kzVQ^{Ro)\t ^ A;r. For axisymmetric modes with negligible vertical 
shear, this is always satisfied. However, for m ^ disturbances, or flows with non-negligible k^VQ^, 
spatially localized wavepackets maintain a characteristic radial wavenumber for only a limited 
time, altering their spatial pattern because of the background shear. For a wavepacket with initial 
wavenumber /sr(0), the radial wavenumber at time t becomes, upon inclusion of —t times the 0{r) 
term in lo in the phase, 

knit) = kn{0) - {mn' + k,v'o,)t. (40) 

Thus, for example, with k^ = 0, the pitch tanp = m/Rkji of a spiral pattern changes by a fraction 
e = \dk-R\/k-R over time t = ekji/\mQ'\. If Rk^^ S> m, the pattern changes slowly compared to 
the orbit time. Among nonaxisymmetric disturbances, the wavepackets with low m/Rku have the 
largest temporal range for which they remain close to normal modes of the system. In the following 
two sections, we present solutions for the growth rates of unstable disturbances determined from a 
local modal analysis (i.e., producing solutions u = I4^(k, i?o)), with the understanding that when 
\mW + kzv'oJ 7^ 0, the modal growth (i.e., ~ e''^'*) with fixed pattern holds only for a limited 
time. In assessing the potential of the non-axisymmetric instabilities we shall identify to affect flow 
dynamics, we will consider their total amplification over times < k^R/m^, for which the spatial 
pattern changes little. 

For simplicity let us define dimensionless variables 

a = u;Ro/va{Ro), = KRo, a^R = ^R-Ro, = RoI^o/vk{Ro), 

n = Rono/vA{Ro), and C = RoVo,{Ro)/va{Ro)- 

Here, C measures the amount of shear in the vertical velocity of the winds, and the "o" subscripts 
in the equilibrium epicyclic and rotation frequencies denote evaluation at the reference point Rq. 
We adopt the power-law equilibria of §2. We now organize the terms in eq. (38) finally to get a 
12th-degree polynomial 

10 

= cj^^ + ^/j(g,z,a;R,Xz,m,Q,C)(T^. (41) 

j=0 

The functional dependences of the coefficient f/s on the parameters are so complicated that it is 
not illuminating to write down the whole expression here. We may obtain more simplified forms 
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for /j 's by sorting out terms and taking the limit of S> 1. Because there are various interesting 
modes which demand different regimes of parameters, however, we keep ah the terms as the general 
local dispersion relation (with xr large) and compute numerical growth rates by solving eq. (41) 
for £7 as a function of the other variables. We present these numerical results in §5. In §6, we 
will classify individual (either unstable or overstable) modes and provide their limiting dispersion 
relations. 

Terquem & Papaloizou (1996) considered only incompressible modes, for disk applications 
where thermal pressure is considerable, by taking a divergence-free displacement vector as a pertur- 
bation eigenvector, thus obtaining a 4th-order polynomial. Our dispersion relation, for applications 
to supersonic MHD winds in which thermal pressure is negligible hence motions are compressible, 
contains information about all possible, either oscillating or unstable, modes of cold MHD winds. 

5. Numerical Solutions of Modal Dispersion Relation 

The derived dispersion relation (41) is a 12th-order polynomial with real coefficients, indicating 
that solutions appear as complex conjugate pairs. Solutions having non-zero real and imaginary 
parts are overstable modes, and solutions with vanishing real parts are unstable modes. In this 
section, we present both types of solutions which are consistent with the local analysis by fixing 
= 10. As wc shall show, both higher values of xr (more spatially localized) and lower values of 
xr (less spatially localized) give qualitatively the same family of solutions as with xr = 10. 

5.1. Axisymmetric Modes of Instabilities 

First, we consider the axisymmetric case with m = 0. For 4 selected sets of parameters, 
wc plot the real and imaginary parts of the unstable and overstable modes in Fig. 4. A Keplerian 
rotation gradient with = is assumed and vertical shear is neglected except in Co. We take Q as 
arbitrary rather than using the relation (6), by allowing that the gravitational force from a central 
object also contributes to the equilibrium rotation velocity. Then, from eq. (5), the normalized 
angular velocity becomes 

= Ac^ = cos^ i - + Gr, (42) 

where Gr = gR{Ro)Ro/v\{Ro) > is the normalized gravitational acceleration. Note that for 
1 + q > (i.e., magnetic fields decreasing outward), equilibrium solutions with i approaching 90° 
require non-zero gravity (because hoop stresses do not confine a primarily-poloidal flow). Also note 
that as Gr strengthens, the initial equilibrium is maintained by the balance between centrifugal 
and gravitational forces, implying that the magnetic force is negligible. 

The behavior of the solutions shown in Fig. 4 (and similar behavior for other parameters) 
allows us to identify 4 different axisymmetric mode families: a toroidal resonance mode (TR), 
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an axisymmetric toroidal buoyancy mode (ATB), a poloidal buoyancy mode (PB), and a Balbus- 
Hawley (BH) mode. One (TR) of these is an overstable mode and the others (ATB, PB, and BH) 
are purely growing modes. 

Fig. 4a and 4b correspond to a disk wind at large distance from the source, where magnetic 
fields are doniinantly toroidal (small i) and centrifugal force balances magnetic force (small Gr), 
while Fig. 4c and 4d correspond to an accretion disk or inner part of a wind where magnetic 
fields are poloidal (large i) and centrifugal force is balanced by the gravity from a central object 
(relatively large Gr). In each frame, solid and dotted lines represent the imaginary and real parts 
of the frequencies, respectively. Fig. 4a shows the TR mode which splits into two branches in the 
presence of (arbitrarily small) poloidal magnetic fields (Fig. 4b and 4c). This TR mode is not 
a generic instability mode because it has a far larger real part (associated with ordinary MHD 
oscillations), indicating an overstability. With the presence of poloidal field components, there 
exist two different types of buoyancy modes, namely ATB and PB modes. When the pitch angle 
of the magnetic field is relatively small, the buoyancy instabilities are driven by the interplay of 
the centrifugal force with the hoop stress of toroidal fields, so we call these ATB modes. Since, 
as explained in §6.1.2, the ATB modes need non-zero poloidal fields as well to be unstable, they 
disappear when i = 0. On the other hand, the PB instability modes arise when the fields are 
predominantly poloidal so that the pressure gradient forces of poloidal fields and the gravity from 
a central object are main driving forces, similar in character to the Parker instability. ATB and 
PB are pure instability modes with Re((7)=0, as shown in Fig. 4. These instability modes operate 
even in the arbitrarily high-fc^ limit because of our cold MHD assumption; otherwise sound waves 
would stabilize short wavelength perturbations, as they do in the Parker instability. Fig. 4d shows 
that the BH instability mode appears when Gr ^ 1, corresponding to dynamically weak magnetic 
fields in the equilibrium; with reduced Gr (also shown in Fig. 4d) , BH is stabilized by radial MHD 
wave motions when xr is large. As discussed in §6.1.3 and §7.1, one interesting finding in our work 
is that the compressible axisymmetric BH mode is strongly suppressed even for Gr large when the 
toroidal field is sufficiently strong; in a cold, Kepler-rotating MHD flow, it is fully stabilized when 
the pitch angle i < 30° (see §6.1.3). 

Fig. 5 shows how the characteristics of unstable/overstable modes change as i and Gr vary 
for the fixed values of = 4, xr = 10, and 9 = C = 0. For a pure toroidal field configuration with 
i = 0°, we observe only overstable TR modes that are almost independent of Gr. As i increases, 
ATB emerges but is stabilized by rotation with Gr large. When i = 45° and q = 0, the buoyancy 
mode disappears because with these parameters the net force from the background magnetic fields 
vanishes (cf. eqs. [5], [6], and [42]). When i > 30° the BH mode strengthens as Gr increases. This is 
because in our normalization higher values of Gr correspond to weaker equilibrium magnetic fields, 
with which the BH instability operates efficiently. At a pure poloidal configuration of magnetic 
fields, BH and PB modes remain unstable (Fig. 5d). Dotted lines at very small Gr in Fig. 5c and 
5d mark the minimum value of Gr, available for given values of q and i, below which no initial 
equilibrium exists (cf. eq. [42]). 
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5.2. Non-Axisymmetric Modes of Instabilities 

When non-axisymmetric perturbations are applied, the cold MHD system responds with 3 
more modes which are either unstable or overstable, in addition to the axisymmetric modes. We 
shall refer to these as non-axisymmetric toroidal buoyancy (NTB), geometric poloidal buoyancy 
(GPB), and poloidal resonance (PR) modes. The PR modes are MHD waves which have non-zero 
azimuthal wavenumbers and become overstable when there is a radial gradient of the axial field, 
analogous to the TR modes. In addition to the above modes, systems with toroidal magnetic 
field configurations and non-zero sound speed are also subject to significant non-axisymmetric 
magnetorotational instability (NMRI) modes. Unlike the other three non-axisymmetric modes, the 
NMRI mode arises due to a differential rotation with dQ /dR < 0, where O is an angular velocity of 
the rotation. More than anything else, the fact that the NMRI in B^-dominated systems needs a 
finite sound speed to be unstable distinguishes it from the axisymmetric Balbus-Hawley instability, 
which can be unstable regardless of temperature for purely axial fields. As we shall discuss later, the 
basic mechanism for the onset of NMRI is quite different from that of axisymmetric BH instability. 
We reserve the discussion of NMRI modes for §7, concentrating here on numerical results for our 
basic cold MHD system. 

Fig. 6 shows the unstable and overstable solutions of the dispersion relations for two combina- 
tions of selected parameters: Fig. 6a corresponds to a disk wind with small i and small Gr, while 
Fig. 6b is for the near-disk case with large i and large Gr. We assume a Keplerian rotation and 
take an arbitrary Q once again using eq. (42). For all cases, we chose = 10, g = = 0, and 
m = 1, and confirmed that changes to these parameters do not appreciably affect the qualitative 
results. Solid and dotted lines in Fig. 6 represent imaginary and real parts of the normalized wave 
frequencies, respectively. When i = Gr = 0, Fig. 6a shows the presence of unstable NTB and over- 
stable TR modes which are split by the non-axisymmetry. NTB modes are nearly like PB modes 
in their physical basis and have an almost constant growth rate over a wide range of Xz- But they 
depend sensitively on the logarithmic gradients of the density and magnetic structures (i.e., q; see 
Fig. 7). For an intermediate value of i, both TR and PR modes coexist. At some wavenumber x^, 
they combine to simply vanish, but overall they give rise to complicated behavior of lm.{a). When 
i = 90°, we observe three unstable GPB, BH, and PB modes, and overstable PR modes (Fig. 6b). 
GPB modes are driven by a buoyancy force together with the geometrical effect. Note that the 
real parts of PR modes are linearly proportional to the vertical wavenumber x^, as TR modes are, 
indicating that they are really overstable modes. Since ,tr S> m, however, there exists only a small 
contribution from non-axisymmetric effects to the axisymmetric BH and PB modes (cf. Fig. 4d). 
Remember that when m ^ xr, the normal mode assumption rapidly breaks down because such 
high-m modes lose their spatial pattern very quickly; we investigate m S> xr cases using different 
methods in §§7 and 8. 

Fig. 7 shows how the characteristics of the buoyancy modes change as x^, i, and q vary. When 
Gr=0, an initial equilibrium exists only for a limited range ofi< icrit = cos~^ •\/(l -|- q)/2 from eq. 
(42), with toroidal field components dominating over poloidal field components. In Fig. 7, therefore. 
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the unstable modes with i < icrit correspond to toroidal buoyancy modes, while poloidal buoyancy 
modes have i > ia-it- Generally speaking, with the assumption of extremely cold medium, smaller- 
scale buoyancy modes with high have larger growth rates. When i is very small, as seen in Fig. 7, 
ATB modes are stable because they need the aid of poloidal fields to be unstable, while NTB modes 
become unstable for all i < icrit • This reflects the physically different driving mechanisms between 
ATB and NTB instabilities. PB modes become more unstable with higher q (steeper background 
gradients), while ATB/NTB modes arc more efficient with smaller q. Greater instability is simply 
associated with higher background magnetic force in the respective cases (cf. the initial equilibrium 
condition [42]). 

The A;R-dependence of the unstable/overstable modes are summarized in Fig. 8. Here we fix 
Xz = 2 for all cases and choose q = 0.8, Gr = 0, and small i for Fig. 8a and 8c, corresponding 
to disk wind-like systems, and q = 0, Gr, = 5, and large i for Fig. 8b and 8d, corresponding 
to accretion disks or disk winds near their sources. The BH instability modes are completely 
suppressed by MHD waves when xr > 3; we will show that this is consistent with the prediction of 
the asymptotic dispersion relation. All the other modes extend with smaller growth rates to larger 
xr, with lm{a) -Tj^^, which we will show agrees well with the asymptotic dispersion relations (43), 
(45), (50), and (52) for the PB, ATB, TR, and NTB modes, respectively. For the PR modes the 
asymptotic dispersion relation (56), showing lm(c7) ^ x-^ , is valid only when Rkz 3> m, which is 
not consistent with the parameters adopted in Fig. 8d. When A;r 3> fcz ~ m/R, one can confirm 
analytically that the PR modes also behave as Im{a) ~ x^^. In the shearing- wavelet point of 
view with cq. (40), Fig. 8 shows that kinematic shear arising from the background flows ultimately 
stabilizes both unstable and overstable modes, as A;r grows secularly increases in time. Although 
the local approximation breaks down if .tr is not large. Fig. 8 indicates that the BH mode exists 
and may show interesting behavior for small xr. In addition. Fig. 8 also suggests larger growth 
rates when xr is small for other modes, although the assumptions of this section of a radially-local, 
slowly-changing pattern are not self-consistent when xr is small. To study dynamical growth of 
disturbances which occurs when xr <C m, we use direct integrations of the shearing-sheet equations. 
We present these results in §8.2 (for the NTB modes) and §§7.2 and 7.3 (for the NMRI modes and 
generalized MRls). 



6. Mode Classification 

The cold MHD system we are investigating has 8 distinct local modes with lm{a) > 0. Some 
of them (TR and PR) have larger Re((7) corresponding to overstability, while the others (PB, ATB, 
BH, NTB, GPB, and NMRI) have negligible Re(cr), indicating pure instability. The NMRI modes 
do not appear in the numerical solutions because of the cold MHD assumption we made. Detailed 
discussion of the NMRI modes will be separately given in §7.2. In this section we describe the 
physical natTire of the individual cold-fluid modes and present the respective dispersion relations 
under some limiting approximations. 
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6.1. Axisymmetric Modes 

6.1.1. Poloidal Buoyancy Mode 

Consider a system with pure axial fields. If gravitational forces are large, then they may balance 
the combined outward radial centrifugal force and pressure gradient force of outward-decreasing 
Boz{R)', otherwise, if = 0, then the strength in the magnetic fields must increase outward for 
an equilibrium to exist. In an initial state, at any point in the system the magnetic pressure 
force acting outward is balanced by the difference between gravity and the centrifugal force acting 
inward. If perturbed, a denser fluid element experiences reduced magnetic forces but unchanged 
centrifugal and gravitational forces per unit mass, and thus it would tend to sink radially inward 
dragging the field line with it; a lighter fluid element would correspondingly tend to float outward. 
Then, in a frozen-in-field condition, the neighboring gas finds itself on sloping lines of force and 
thus slides inward to add its weight and to cause field lines to bend more, expediting the instability. 
This poloidal buoyancy mode is analogous to the Parker instability (Parker 1966), with the driving 
force role of external gravity in Parker's instability replaced by combination of gravity and the 
centrifugal force in the PB. The PB mode can occur for both axisymmetric and non-axisymmetric 
disturbances. 

Putting BQff, = m = and considering short wavelength perturbations with vj^^k^ S> in eq. 
(38), one can find the dispersion relation for the poloidal buoyancy mode is 

_ o ^ _ fl + q\ ' vlM ^ _ f \9K-m'n kl 

\ 2 ; R\kl + kl) \ < Jkl + kl ^^'^ 

Eq. (43) states that there is no preferred length scale as long large. However, the 

inclusion of thermal effects would stabilize the PB mode with shorter wavelengths, as in the Parker 
instability^ Also, sufficiently large /cr /cz stabilizes the mode. 



^Also by taking a local approximation and by neglecting density stratification and the eflfects of thermal and 
cosmic ray pressures, one can simplify eq. (III. 12) of Parker (1966) to get the asymptotic (k oo) dispersion relation 



where g is the gravity perpendicular to the galactic plane, iiah is the Alfven speed of initial fields parallel to the 
galactic plane, and fcji and k± are perturbation wavenumbers in the respective directions parallel and perpendicular 
to the galactic disk and magnetic field. Comparing the above with eq. (43), we may write ^efi = ffR — for the 
PB modes, with wavenumber correspondence fez <-> feii and fen <-> fex. 
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6.1.2. Axisymmetric Toroidal Buoyancy Mode 

Now consider a system with weak poloidal but strong toroidal field components and negligible 
gravity. When magnetic fields are predominantly toroidal (i.e., when i < cos~^ ^/(T+q)/2 from 
eq. [42]), an initial equilibrium state is maintained by the balance mainly between the centrifu- 
gal force acting outward, and magnetic hoop stresses which act inward. With sinusoidal density 
perturbations with kz imposed on the equilibrium, a heavier blob of material would tend to float 
radially outward under the action of unchanged centrifugal forces per unit mass but reduced speciflc 
magnetic forces; a lighter element would correspondingly tend to sink. The radial motions of the 
heavier and lighter blobs are in opposite directions and thus cause the poloidal field lines to bend, 
creating radial perturbed fields. 

The azimuthal fiuid motion is slightly accelerated by the tension force exerted by the initial 
toroidal and the perturbed radial fields (cf. Bq^Biji/R term in eq. [28], associated with spiral 
magnetic field line projections in the 2;=constant plane). This causes the initial poloidal component 
of field lines to bend now in the azimuthal direction, creating bands of perturbed azimuthal fields 
with signs alternating in the z direction. The resulting total azimuthal fields arc distributed in such 
a way that the heavier (lighter) blob in the initial perturbation has a lower (higher) toroidal field 
strength. Induced motions due to the vertical magnetic pressure gradient force carry the matter 
from under dense to over dense regions, closing the loop and amplifying the initial perturbation. 

By setting m = vqz = and taking the v\k'^ S> limit, we obtain from eq. (38) the following 
dispersion relation 



= 0^- 



2 -"0z 
^0 



oj^ - 4nvAct>VAzFBkzU; - vj^vX^F^k^, 



(44) 



with Fb = (cos i — (1 + q)/2)/R. When vx(t, = (i.e., with pure poloidal fields), eq. (44) 
immediately recovers eq. (43), the limiting dispersion relation for PB modes. On the other hand, if 
vat. = 0, there is no unstable ATB mode, clearly demonstrating that ATB modes operate by bending 
poloidal field lines. Prom Fig. 5b, we note that ATB modes are stabilized by rotation (larger Gr 
corresponds to stronger rotation). It can be shown from eq. (44) that when Rk^, Rk^ ^ u), the 
critical wavcnumbers are (fc^ + /cR^)crit = — {d^'^/d In R) /{v\ svo? i) , below which the system is stable 
against the ATB modes. For k^, 3> k^^cnt-, eq. (44) is further reduced to 
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(45) 



Since ATB instabilities are axisymmetric modes, they can persist without being disturbed by the 
kinematic growth of A;r due to shear, if v'q^ = 0. 

Among well-known plasma modes, the pinch or sausage mode of a plasma column is most 
similar to the ATB in overall geometry and effect. Both axe axisymmetric and require the radial 
tension force from predominantly toroidal magnetic fields to drive the instability. For both the 
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plasma pinch mode and the ATB of cold cylindrical winds, the net effect is that matter tends to 

be ejected radially in bands alternating with contracting magnetic field loops. However, in pinch 
modes, the plasma is generally unmagnetized and surrounded by external toroidal fields, and axial 
fields tend to suppress the instability. In the ATB, on the other hand, internal toroidal magnetic 
fields permeate the fluid, and non-zero axial fields are required for instability. 



In the presence of axial magnetic fields, a differentially rotating disk is unstable to an ax- 
isymmetric incompressible perturbation (Balbus & Hawley 1991; see also Velikhov 1959 and Chan- 
drasekhar 1960). Because this Balbus-Hawley instability^ has a rapid growth time (comparable to 
the local rate of rotation) and exists for arbitrarily weak magnetic field strength, it is believed to 
provide a powerful mechanism for the generation of the effective viscosity in astrophysical accretion 
disks. Through numerical simulations, Hawley &: Balbus (1991) argued that the roles of compress- 
ibility and toroidal fields are not significant as long as the total field strength is subthermal. Also, 
Blaes & Balbus (1994) studied the effect of toroidal fields on the compressible axisymmetric BH 
instability and showed that toroidal fields do not modify the instability criterion, while reducing 
growth rates slightly if vji^ff, < Cg. We find the striking result that under extremely cold conditions 
(i.e., VA ^ Cs), compressibility prohibits the axisymmetric BH instability from occurring if the 
toroidal fields arc as strong as the poloidal fields. 

By taking the weak magnetic field limit (0, ^ vp^kji, vp^kz) and m = 0, we obtain from eq. (38) 
the following dispersion relation for the compressible axisymmetric BH instability in a cold MHD 
flow 



With an n oc rotation profile, eq. (47) implies that if sin^z < 1 — a/2, we anticipate no BH 
instability in a cold flow. For a Keplerian rotation law with a = 3/2, for instance, no axisjnnmetric 
BH instability occurs if i < 30°!; when the magnetic field strength is superthermal, the inclusion of 
toroidal fields tends to suppress the growth of the BH instability. With a steeper rotation profile 
(as would occur, for example, in winds from boundary layers), there is an increase in the range of 
i for which a system is BH-unstable. 

We defer the full discussion on the BH instability until §7.1, where we explicitly include pressure 
terms in the dynamical equations. 



6.1.3. Compressible Balbus-Hawley Mode 




(46) 



4z(fcz + 4) + - 4^^^ sin^ i < 0. 



(47) 



^Often referred to magnetorotational instability, or briefly, MRI. 
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6.1.4- Toroidal Resonance Mode 

Consider a system having pure toroidal fields without rotation. If the initial fields are homo- 
geneous in space, magnetosonic waves, driven solely by magnetic pressure (with the assumption of 
the cold medium), would propagate without any interruption in the plane whose normal is per- 
pendicular to the magnetic field direction. In an inhomogeneous medium, however, MHD waves 
no longer maintain a sinusoidal planform, and the characteristics of the waves change through the 
interaction with the background medium. The amplitudes of the waves may sometimes increase as 
they propagate, or sometimes they may become evanescent and decay at a resonance, or even may 
be trapped between two resonance points (cf. Rae & Roberts 1982). In such a strongly structured 
medium, the classification of MHD waves is not in general possible. 

Our local treatment of MHD waves can provide some insight on the amplification or evanescence 
of propagating MHD waves in a structured medium. For Up ^ and Bqz = $7 = 0, the local wave 
equation (35) can be simplified as 

d^^R dlnw|(i^R w| 

iB^ - -iR-ltR + = ^^^^ 

Again we define = (a;|)^/^*, then eq. (48) takes the form of eq. (16) with ^r = K{R) defined 

by 

(49) 

where vxcj) oc R~^/'^ was assumed. Thus, considering limiting cases of and A;r, one can find the 
dispersion relation for this mode near the resonance frequencies (i.e., a) k, fA(/>fcz) 




' v\^kl [1 + (3/4)V2e±«/2(ii/cj^)-i] , for Rk^ » Rk, » 1, 

(50) 

v\^kl [1 + (3/4)V3e±2-i/3(^fc^)-2/3j ^ for ^ ^ 1^ 

showing that the imaginary part of toroidal resonance mode vanishes quickly as A;r — > oo, while 
its real part gets bigger as k^ increases. Therefore, it is not adequate to regard TR modes as a 
true local instability mode. Though the TR mode is not a local instability, it suggests potential for 
waves to have global instabilities in which the magnetosonic resonance (a) = VA(f>kz) plays a similar 
role to the Lindblad resonance in rotating disks. Thus, waves of fixed frequency propagating with 
a radial component of k into their magnetosonic resonances may be amplified or reflected. 

The modification of traveling waves due to the inhomogeneity of the medium is mediated 
through the magnetic pressure. A similar effect would occur when hydrodynamic waves propagate 
into an inhomogeneous medium. 
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6.2. Non-Axisymmetric Modes 

6.2.1. Non-Axisymmetric Toroidal Buoyancy Mode 

The non-axisymmetric toroidal buoyancy mode is very similar to the PB mode in its physical 
mechanism, in spite of the different field geometry. For toroidal-field dominated cases, an equilib- 
rium can exist with the net inward magnetic stresses balancing the outward centrifugal force. When 
the system is perturbed non-axisymmetrically, the instability would develop similarly to PB modes, 
as described in §6.1.1. For the B^-dominated case, however, over-dense regions float outward and 
under-dense regions sink, because the inward magnetic forces are enhanced when the density drops. 

Setting Bqz = ^Oz = and assuming Rkz <C m, Rkn, one can show that the general dispersion 
relation (38) is reduced to the following quartic equation in terms of u 



= iu^- 



vl, + 4) + - 4^^^FB.i,^ - vi,Fi , (51) 



with Fb = (1 — q)/2R, The negative last term in eq. (51) guarantees the existence of unstable NTB 
modes. The third term (caused by the coupling of the rotation with the background fields) tends 
to stabilize NTB modes. Thus, if 

there is no unstable NTB mode. Note that for wind equilibria with i?oz = 0, from eq. (6) the 
RHS of the above equals 3(1 — q)/2R^\ thus the NTB instability will be present at all m when 
1/3 < ^ < 1. 

In the limit of large m, we obtain the asymptotic dispersion relation for the NTB mode 

V 2 y R-^im^ + R^kl) \ vl^ Jm^ + R^kl' ^ ' 

Here for the second equality, eq. (6) with i = 0° is used. Eq. (52) is akin to eq. (43), the dispersion 
relation for the PB mode, and to eq. (45), the dispersion relation for the ATB mode, reflecting 
the common origin in buoyancy forces of all three. In fact, Fig. 7 clearly shows how the various 
buoyancy modes extend and smoothly join at intermediate pitch angles. 



6.2.2. Geometric Poloidal Buoyancy Mode 

Now suppose a system with pure vertical fields. When perturbed azimuthally, a fluid element 
becomes over dense and tends to move inward due to the decreased background magnetic pressure 
force per unit mass if < g < 1. This geometrically converging motion of fluid increases density 
and field strength by factors of (1 — q) and (1 — q')/2, respectively. On the other hand, the magnetic 
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field enhancement induces diverging motions of the fluid in the azimuthal direction by building up 

a pressure gradient, tending to lower the density. The net effect of these two processes is a density 
increase by a factor of (1 — g)/2, accelerating the inward motion of the heavier element. When 
m 3> Rkz, the dispersion relation for this GPB mode is found to be 



When q = 1, there is no instability. This is because the initial configuration of the density and the 
field is such that the mass and magnetic flux contained in a thin ring with the thickness dR and 
the radius R are constant over R, and no gain from the geometrical effect is possible. 



The physical basis for the poloidal resonance mode is quite similar to that of the toroidal 
resonance mode. The only difference between them is the background field geometry. In the 
presence of pure axial fields, MHD waves with non-zero m are easily influenced by radial magnetic 
pressure gradients. 

To derive the dispersion relation near the resonance frequencies (i.e., up ^ 0), let us suppose 
a system with pure axial flelds and neglect the vertical velocity shear. The system is also assumed 
to rotate slowly enough that the effects of rotation may not be important in the wave dynamics 
(i.e., mvAz ^ RO,). For —>■ v\^{k^ + m^/R^), we are left from eq. (35) with 



We now define = (a)p/a)^)^/^* to simplify eq. (54) into eq. (16) with kn = K{R) defined by 




(53) 



6.2.3. Poloidal Resonance Mode 




(54) 




(55) 



where we took the limit of Rk^, 3> m and assumed vaz R Solving eq. (55) for two limits of 
/cr, we obtain the dispersion relation near the resonance frequencies 



' v\^kl [1 + e^^/^{ma^^fl^{R^kKK)-'^'^] , for Rk^ » RK » m, 
up' = I (56) 
^ [1 ± e-^'/^{mas^)^/^{Rkz)-^] , for Rk^ » Rk^ » m. 



where a^h = 1 ± 3mQ/vAzkz, showing again a rapidly declining imaginary part as /cr increases, at 
which the local approximation is valid. Thus, just like TR modes, PR modes are not strictly local 
instability modes. 
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7. Magnetorotational Instability 



7.1. 



Axisymmetric BH Instability 



In the preceding section, we briefly discussed the axisymmetric BH instability in a cold, dif- 
ferentially rotating medium and found that the BH instability can be suppressed by the azimuthal 
component of magnetic fields, if the medium is cold enough. Incompressibility has generally been 
adopted in the study of the BH instability in an accretion disk on the grounds that in such a 
system the magnetic fields are subthermal and thus acoustic waves can maintain the incompress- 
ible condition over many rotation periods. For magnetocentrifugally driven winds, however, sound 
waves play a minor role in controlling the dynamics and thus the incompressible approximation 
is inapplicable. In addition, since an initial equilibrium is attained through the balance between 
the centrifugal and magnetic forces (cf. eq. [5]), the Alfven crossing time scale is comparable to 
the rotation time scale (^a ~ in this case, the fields are not weak and the unstable range of 

wavenumbers becomes narrow. 

We generalize the previous discussion of the axisymmetric compressible BH instability by 
explicitly including the thermal pressure terms in the momentum equation and exploring the role 
of compressibility to the development of the Balbus-Hawley instability. We consider a cylindrical 
flow threaded by both vertical and azimuthal magnetic fields, ignoring the radial variations in the 
initial configuration except 0, = Q{R) and neglecting the vertical velocity. We assume the medium 
is isothermal and take the WKB (Rkz ^ 1) approximation. Through the standard approach to 
linear analyses, we arrive at the dispersion relation for the compressible version of the BH instability 



where Cg is the isothermal sound speed of the medium, u!\ = — v\^k^, and /(o;^) is defined by 



Eq. (57) is a sixth-order equation for to with only even terms. When kji = 0, eq. (57) is identical 
to eq. (64) of Blaes & Balbus (1994) or eq. (99) of Balbus & Hawley (1998). Now let us take the 
two opposite limits of Cg to obtain the following dispersion relations 



(57) 



ui-{n^ + vlkl + vl^kl)ujl + {K't,l^-A^''vl,)kl = 0, for ^ 0, 
{l + kl/kDuji-K'ul-A^^vl^kl = 0, for Cg^oo, 

and the corresponding instability criteria"^ 



(58a) 
(58b) 



^Az(^z + ^r) + - 4r2^ sin^ i < 0, for Cg 0, 
^'AziM + 4) + dn^/dlnR < 0, for oo. 



(59a) 
(59b) 



*In fact, from eq. (57) the formal instability criterion (59b) is generic for any value of Cs 0; it may be written 
as v\^{k'i + k%) + — 4f2^ < 0. However, when Cs/va ^ 1, growth rates for small i are very low. 
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Note that eq. (58b) is the same as the original dispersion relation of the incompressible BH 
instability (eq. [2.9] of Balbus & Hawley 1991 without the Brunt- Viiisala frequency). The instability 
criterion (59a) in the extremely compressible limit depends explicitly on the local pitch angle, 
showing that as i departs from 90° the instability becomes gradually confined to smaller values of 
kz- For a cold Keplerian flow, no instability occurs when i is smaller than 30°. 

To examine what role thermal pressure plays to the growth of the BH instability and why 
the instability criterion depends on i, we plot the unstable solutions of eq. (57) as functions of 
gA = (k ■ va)/^^ (= kzVAz/^ for the axisymmetric case) and l3 = / v\ in Fig. 9. For the time 
being, we confine our discussion to the A;r = case. When i = 90°, the instability criterion 
from eqs. (59a,b) is v\k'^ < —dn'^/dlnR and the growth rate is independent of P, implying that 
the compressibility does not alter the instability (Fig. 9a). This can be understood as follows: 
when magnetic fields are mainly axial, sound waves propagating along a vertical direction decouple 
completely from the magnetic fields and are undisturbed by rotation. But transverse MHD waves 
which are intrinsically incompressible are influenced by rotation to become ultimately unstable for 
a range of kz when d^/dR < 0. Therefore, i = 90° is a very special case. 

On the other hand, when both vertical and azimuthal fields are present, toroidal perturbed 
fields generated by an initial azimuthal displacement or by sheared motion following a radial pertur- 
bation of the initial axial fields tend to cause vertical oscillations, but in a cold assumption, mainly 
due to the magnetic field gradient terms, —BQ^{dBi(f,/dz). This oscillatory vertical motion tries 
to distribute the perturbed fields as uniformly as possible, thereby tending to suppress the growth 
of the disturbances. However, the vertical magnetic pressure gradients are not strong enough to 
create significant vertical motions if thermal pressure is large: a compressed region tends to expand 
vertically but with little change in the strength of the toroidal fields, thus providing a favorable 
condition for the development of the BH instability. This explains why higher /3 cases have higher 
growth rates at fixed i, and why the growth rate decreases as i decreases at fixed /3 (Fig. 9b~9d). 

Although the instability criterion (59b) is completely independent of the strength of the az- 
imuthal fields provided that /? 7^ 0, indicating as noted by Blacs &: Balbus (1994) that to all orders, 
azimuthal fields do not modify the stability criterion, the corresponding growth rates drop progres- 
sively as (3 decreases Hi ^ 90°. When /3 > 1, any change of an inclination angle i from 90° does not 
bring significant reductions in growth rates, implying that the characteristics of the instability are 
essentially the same as the pure poloidal case. If /3 ^ 1, however, we observe dramatic stabilizing 
effects from toroidal fields, as illustrated in Fig. 9. 

A few comments should be devoted to the effect of A;r. Radial wave motions do nothing but add 
another restoring force to perturbations. This in turn means that thermal pressure has a stabilizing 
influence on the growth of the BH instability. Thus there are two competing processes of thermal 
pressure: thermal pressure associated with vertical wave motion promotes the BH instability, while 
thermal pressure controlling radial motion opposes it. It turns out that for i 7^ 90° the former 
process always dominates. For i = 90°, only the latter effect exists, giving higher growth rates for 
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smaller /3, when A;r 7^ 0. 

Notice the stabihzing effect of A;r in eqs. (59). If the background vertical flow has significant 
shear, the local radial wavenumber would increase with time (cf. eq. [40]), suppressing the instability. 
Thus when 7^ 0, the compressible BH instability will exhibit a transient growth, as must happen 
to all modes if kzv'o^ / and/or m^l' / 0. 

In conclusion, we have found that compressibility has a stabilizing effect on the axisymmetric 
BH instability. Even though its effect is small if the sound speed is super- Alfvenic, compressibility 
must be considered whenever the Alfven speed is comparable to or even exceeds the thermal sound 
speed, as is expected in winds and also in disk coronae (cf. Miller & Stone 1999). 

The above discussion applies only for axisymmetric perturbations. It was also found that 
an accretion disk with purely toroidal fields is subject to non-axisymmetric instability (Balbus & 
Hawley 1992; Terquem Sz Papaloizou 1996), but we will show in the following section that the 
physical role of compressibility in that case is completely different, in spite of the same quantitative 
instability criteria. 

7.2. Non- Axisymmetric MRI: Coherent Wavelet Analysis 

Balbus & Hawley (1992) found that a differentially rotating disk of incompressible fluid with 
embedded toroidal magnetic field is unstable to non-axisymmetric perturbations. Adopting shearing 
sheet coordinates (sec below), they integrated a set of the perturbed equations and showed that 
perturbations with an intermediate azimuthal wavenumber m can exhibit transient, but enormous 
growth over a time scale of several percent of 

An alternative approach was taken by Terquem & Papaloizou (1996) to study a similar insta- 
bility to that identified by Balbus & Hawley (1992). They solved the problem using the local WKB 
approximation. They started from a general compressible equation of state, but subsequently they 
supposed divergence-free poloidal Lagrangian displacements, which made their treatment essen- 
tially incompressible. They derived a sufficient condition for the instability which is exactly the 
same form as that of axisymmetric BH instability (i.e., dQ'^ / din R < 0). Noting that azimuthal 
shear is the main driving mechanism and bending of the field lines provides a stabilizing restoring 
force, they suggested the non-axisymmetric instability of toroidal magnetic fields might resemble 
the original BH instability. 

We argue in this work that the underlying physical mechanisms for non-axisymmetric toroidal- 
B MRI (which we refer to as "NMRI") and axisymmetric poloidal-B MRI (which we refer to as 
"BH") are in fact quite different from each other. In this section, we analyze the NMRI by looking 
at "coherent wavelet" solutions in which every physical variable, localized in both space and time, 
oscillates or grows with the same space-time dependence, and provide quantitative results in detail. 
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1.2.1. Localization in Space and Time 

We begin by considering a shearing, rotating disk with uniform density, and magnetic fields 
with only an azimuthal component. We ignore any unperturbed vertical motion in the medium. 
We include thermal pressure effects with an isothermal equation of state to obtain the explicit 
dependence of the NMRI on the tempcratTire, but neglect effects of cylindrical geometry. This 
configuration is the same as Balbus & Hawlcy's (1992), except that they considered only the in- 
compressible case with the Boussinesq approximation, and allowed for vertical equilibrium gradients 
yielding buoyant oscillations. Adopting the shearing sheet coordinates {R, (p, z) such that R = R, 
^ = (f) — Q,[R){t — to), and z = z (Goldreich &; Lynden-Bell 1965; Julian & Toomre 1966; Balbus 
&; Hawley 1992), we consider the time development of an initial plane- wave disturbance which 
preserves sinusoidal variation in the local rest frame of the equilibrium shearing, rotating flow 

XI {R, 0, z, t) = XI {t)e''^4>+ik.z+ikn(to)R^ (60) 

where k^{to) is a radial wavenumbcr at a fiducial time t = to. The linearized form of the MHD 
eqs. (1)^(4), can be written in dimensionless form as 

^ = -QRUm - QmUlcf, - QzUu, (61) 

^^^^ = 2ui0 - QmbR + qniPa + b^), (62) 



dr 



du 



2 



dr 2^2 
duiz _ 
dr 

dbR 

dr 



uiR + Pq„ia, (63) 
qrabz + QziPce + b^), (64) 

QmUlR, (65) 



db^ dlnn 

-dF = dl^^^ - ^^^^^ - ^^^^ 

— = qraUl2., (67) 
dT 

where the dimensionless Lagrangian derivative is denoted by 

d^-nm^W ^^^^ 

In eqs. (61)~(68), all perturbed variables are dimensionless and defined by a = pi/po, Ui = 
i'Vi/vAtj,, b = Bi/i3o0, and r = tfl, and dimensionless parameters are /? = Cs/v\^, q^ = VA<t)m/Rn, 
qz = VA(t>kz/^, and 



, . VAct>kR{t) VA4, 



kR{to) - m(t - to)-^ 



.mt^—--Ta ^ (69) 
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where the third equaUty holds when to = —k^{to)lm^'. Eqs. (65)~(67) yield the divergence free 
condition for the perturbed magnetic fields. 

Since has a r-dependence, the linear system of eqs. (61)~(67) does not form an eigenvalue 
problem; kinematics of shear wrap a given disturbance by increasing its radial wavenumber linearly 
with time. In the original shearing sheet formalism, the fate of a system exposed to perturbations 
is analyzed through direct integrations of linearized equations over time. In doing so, one may 
observe transient amplification or decay of applied disturbances depending on their stability. One 
can say that a system is unstable if some physical variables grow sufficiently over certain time scales. 
The efficiency of instability for a system is identified by computing the response of the system to 
variation of parameters input to temporal integrations. This approach was adopted by Balbus &: 
Hawley (1992) in their identification of the NMRI. 

Here, we instead analyze the NMRI by proceeding one more step from the original shearing 
sheet formalism to find solutions which are localized in time as well as in space. First, we note that 
there exist two distinct time scales: the growth time of instabilities determined by the inverse of 
the dimensionless instantaneous growth rate, 

and the dimensionless shearing time, [dhiq-^/ dT)~ = gR|(/m 

dlnJ)/dlni?|-i typical time scale 
of the linear growth of the radial wavenumber. If the shearing time is much longer than the growth 
time, that is, if (^mlfiln il/dln i?|/(7gR) ^ 1 (the "weak shear limit"), the time dependence of 
in eq. (69) can be neglected, and thus normal mode solutions having an exponential or oscillatory 
behavior can be sought. Shu (1974) applied this technique to investigate the effects of a differential 
rotation on the Parker instability. Also, Ryu &; Goodman (1992) obtained an algebraic dispersion 
relation for the convective instability in differentially rotating disks, by assuming that q^ is time- 
independent. 

Because the convective and the Parker instabilities arise from hydrodynamic and magnetic 
buoyancy effects, respectively, independent of the rotation of a disk, one can always find a regime in 
which the weak shear limit is applicable. In some eases, however, as for example in the axisymmetric 
poloidal or the non-axisymmetric toroidal MRIs with weak magnetic fields, the instabilities result 
directly from a differential rotation with < 0. In such cases, peak growth rates are of the same 
order as rotational frequencies (Balbus & Hawley 1998), and thus the weak shear is not a good 
approximation for these non-axisymmetric instabilities. 

However, we can still look for coherent solutions in which all perturbed variables vary as e^'^ 
with time, provided the variation of the instantaneous growth rate 7(r) over the growth time 'y~^ 
is relatively small, i.e., 

(iln7(T) 



«7(r). (71) 

We refer to the solutions under this approximation as "coherent wavelet solutions" because all phys- 
ical quantities localized in both space and time grow at the same instantaneous rate. If the condition 
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(71) holds, the changes in j{t) can be neglected over a short time interval, and the set of dynamical 
equations (61)~(67) constitutes an cigcnsystcm instantaneously. This is equivalent to the WKB 
method in the time dimension. Since ^^^dln^{T)/dT = — q^iid In Q/ din R){d In din qYi}/{'jqYi), 
eq. (71) is satisfied if either qj^\dlni^/dlnR\/{-yq^_) <C 1 (the weak shear limit), or [dln'y /dlnq^il <C 
1 (instantaneous growth rates are relatively insensitive to the radial wavenumber); the condition 
(71) is less restrictive and in fact is the generalization of the weak shear limit. Of course, we need to 
check the self-consistency of this coherent wavelet approximation by examining a posteriori whether 
resulting solutions satisfy the condition (71). For incompressible media, Balbus &: Hawley (1992) 
mapped the regime of instability in [(k • va)^, |k|/fcz] space using WKB methods similar to those 
we adopt. 



7.2.2. Coherent Wavelet Dispersion Relation 

Upon substituting eq. (70) into eqs. (61)~(67) and applying the approximation (71) so that 
dxi/dr — > 7Xi, one can form a matrix equation 7Q = M.Q, where Q = (a, Uir, ui^, 1x12,6^, b^, 6^)^ 
is a column vector and At is a 7x7 matrix whose components are determined by the coefficients 
of Q in the right hand sides of eqs. (61)~(67). By solving the condition dct(Al — 7X) = 0, where 
I is the identity matrix, we obtain a seventh order polynomial in 7. As a further approximation, 
however, if at least one of the conditions, q^ ^ qr, 7 S> t, or 7r 1, is satisfied, all even order 
terms that depend linearly on r and q^ but are independent of qz, can be neglected compared to 
the remaining terms. The first two conditions apply when the radial wavenumber is not significant, 
either because disturbances are highly localized in the vertical direction {qz S> qr) or simply because 
we are looking at modal behaviors at the time r ~ 0, while the third condition holds when net 
amplification of perturbations is large. This simplification yields a trivial solution 7 = (this 
arises from the fact that perturbed magnetic fields, 5r, 6^, and b^, are linearly dependent via the 
divergence-free condition) and a third-order polynomial in 7^ which is the resulting instantaneous 
dispersion relation for NMRI 



= 7^ + 7^ 



+ Y 



{l+(3)q^+ql + — 



{l + 2P)q'qi + —{ql + {l + P)q^,) 



2 2 I 



dlnQ^ 
d\nR' 



(72) 



where the amplitude of the total wavenumber defined by q'^ir) = qRir)'^ + Qm + 1z is a function of 
r through eq. (69). Combining eqs. (69) and (72), one can evaluate a local, instantaneous growth 
rate at a given time. 

With vanishing magnetic fields and thermal pressure, we would obtain from eq. (72) stable 
epicyclic oscillations. In the limit of strong magnetic fields and no rotation, eq. (72) is immediately 
reduced to (7^+Qm)(7^+(l+/^)9^7^+/^Q'm^'^) = 0) the usual dispersion relations for the Alfven waves 
and the fast and slow MHD waves in a medium embedded with toroidal magnetic fields. In the 
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presence of rotation with non-vanishing but weak fields, however, these Alfven and MHD modes are 
coupled to exhibit generally complex modal behaviors. They can be stable or unstable depending 
on the parameters, but it is always a slow MHD wave that becomes unstable because it has the 
lowest frequency so that there is a plenty of time during which destabilizing forces (centrifugal 
forces for NMRI) act on it. We have instantaneously growing solutions with real positive values of 
7 provided that the last term in eq. (72) is negative. Thus, when fiq^a 7^ 0, the local, instantaneous 
instability criterion in terms of the dimensionless variables for the NMRI can be written 



demonstrating that dTt^ jdhiR = — AfP < is indeed a sufficient condition for the instability 
to arise when the magnetic field strength is negligible (Balbus &; Hawley 1992; Foglizzo & Tagger 
1995; Terquem & Papaloizou 1996). Eq. (73) recovers the results of Balbus Sz Hawley (1992) 
for the instability regime for toroidal-field NMRI. The result of eq. (73) can be compared to the 
poloidal- field BH instability criterion of eq. (59b). Notice that the NMRI (with Bq^ = 0) vanishes 
completely as /3 ^ 0, while the axisymmetric BH instability still exists even when (3 = (see §7.1). 

Although they share the same instability criterion, the operating mechanisms for the NMRI 
of toroidal fields are quite different from the axisymmetric BH instability of poloidal fields. Both 
arise via destabilization of the slow mode. The NMRI mode, just as the poloidal BH mode, 
depends on shear to generate azimuthal fields from radial perturbations of the background fields. 
But there is more to the story. The key mechanism for the NMRI instability lies in the vertical 
MHD wave motions driven by the gradient of the total (initial plus perturbed) azimuthal fields, as 
schematically illustrated in Fig. 10. Since we suppose perturbations which are sinusoidal in both 
vertical and azimuthal directions, the perturbed azimuthal fields are also periodic in both directions. 
Rapid vertical motions with high kz, generated by —Bo^{dBi^/dz) stress, would produce over- and 
under-dense regions which regularly alternate along the azimuthal direction (Fig. 10a). And then, 
azimuthal fluid motions are induced, according to the equation of continuity; from over-dense regions 
to under-dense regions (Fig. 10b). Depending on the direction (t*^) of these induced motions, the 
coriolis and/or centrifugal force would alter the paths, radially inward or outward (Fig. 10c). Under 
the condition of field freezing, these radial motions would produce radial magnetic fields with a 
small amplitude from the background toroidal fields (Fig. lOd). These radial fields would in turn 
be sheared out to generate (positive or negative) perturbed azimuthal fields, due to the differential 
rotation of the background flows. When dQ/dR < 0, the resulting azimuthal fields from initial and 
perturbed ones would be distributed (Fig. lOe) such that they reinforce the applied initial vertical 
perturbations (Figs. 10a and lOf), implying the MRI; the entire system would just oscillate with 
rotation- modified MHD frequencies if dQ/dR > 0. This explains how the NMRI operates. 

When kz is large, the stabilization of the NMRI occurs when a magnetic tension from radially 
bent field lines exceeds the centrifugal or coriolis force (Figs. 10b and 10c). Shear Alfven waves 
with radial polarization can suppress the instability if the field lines are sufficiently strong or if the 
azimuthal wavenumber is large enough, as expressed by the dimensionless parameter outside 
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the parentheses in eq. (73). When -^i qr, on the other hand, MHD waves propagating along the 
radial direction stabilize the NMRI, as indicated by the terms inside parentheses in eq. (73); qR{T) 
clearly reflects the stabilizing effect of the background shear. 



The maximum instantaneous growth rate is achieved when qr 
implies instability if 

1 
2 



0. In this case, eq. (73) 
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d\nR 



Note that gm,crit ~^ -dlnQ.'^ /dlnR, for > 1, while ^m.crit ~^ <l2.y/-dln Q.'^/dlnR, for < 1. 
Numerical solutions of eq. (72) with = are presented in Fig. 11. As both eq. (73) and Fig. 11 
show, the maximum growth rates for toroidal-field background states are attained when — ^ oo, 
which is a sharp contrast to the axisymmetric poloidal-field BH instability that has fastest growing 
mode at moderate ^z's (cf. Fig. 9 and see also discussion in Balbus &; Hawley 1998). But both forms 
of the MRI have the same maximum growth rates at the same = Fig- Hb shows how 

growth rates depend on the sound speed. As the sound speed increases, the medium becomes more 
unstable. This reflects the incompressible nature of the NMRI. Even though the marginally critical 
wavenumber is independent of temperature (for (3^0), the virulence of the instability is greatly 
inhibited as (3 decreases. For q^ ^ 1, one can find from eq. (72) the temperature dependence of 
the limiting growth rate 
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k2(1 + /3) dlnR' 



(74) 



or 7 = qm\/3(3/ (1 + (3) for a Kcplcrian rotation. Eq. (74) gives slopes of the growth rates for small 
qia (Fig. 11). For magneto centrifugally driven winds which are as cold as /3 < 0.01, the NMRI is 
not expected to play a significant role; the growth rate in dimensional units is \/3csm/-R. 

When the medium is incompressible (/3 oo), cq. (72) allows the analytic expression for the 
instantaneous growth rate for pure toroidal- field background states, 
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dlnR 



<0, 



(75) 



otherwise, 



When ^ 1, one can derive the maximum growth rate 7max = \dln^}/d\nR\/2, which is achieved 
when (?m,max = —{dlii^}'^/dliiR)/2 — 7max- It can be shown from eq. (72) or (75) that dj'^/dr ~ 
QmlK/o^ — as — oo. This proves that the coherent wavelet approximation is self-consistent 
for the NMRI with high q^. 



7.2.3. Comparison With the Shearing Sheet Formalism 



In order to compare the coherent wavelet solutions for toroidal-field NMRI with the results 
from the shearing sheet approximation, we directly integrate eqs. (61)~(67) over time, with given 



sets of initial conditions. In Fig. 12, we display the time evolution of all perturbed variables for 
Qm = 0.1, Qz = 1, and /3 = 100, which are the same parameters as chosen for Fig. 3 of Balbus 
& Hawley (1992). We adopt a Keplerian rotation profile in what follows. The initial amplitudes 
are 0.1 for every variable except 6r = 0.01 and bz = 0.4, and the initial f, where the orbit 
number f = t I2it = tQ./2TT, is allowed to be determined from the divergence free condition of the 
initial, perturbed magnetic fields. When f < —20, the system responds with MHD wave motions 
before they start to grow. During this relaxation stage, fast MHD modes having large I^rI are 
nearly longitudinal acoustic waves, affecting wir and a, while ui^ and h^f, are mostly influenced by 
transverse slow modes. As time increases, Ig^l gradually decreases, permitting rotational shear to 
affect the overall dynamics. Once the condition (73) is satisfied, shear drives the slow modes to be 
unstable, following the process illustrated in Fig. 10. Even though the growth of disturbances shows 
a transient nature due to the kinematic growth of qr, the net amplification is about 9 orders of 
magnitude, a bit higher than Balbus & Hawley's result. This is because the integration interval in 
Balbus & Hawley covered a slightly smaller part of the unstable time range. At later time when 
has a large value, the system exhibits stable oscillations with the slow MHD wave frequency. Fig. 
12 also shows the predicted amplification magnitude (thick solid line) from the coherent wavelet 
approximation (see below). 

In Fig. 13, we plot the numerical growth rates for each variable calculated from Fig. 12 based 
on the direct numerical integrations in the shearing sheet formalism, together with the growth rate 
of the corresponding coherent wavelet solution. Here, a dimensionless instantaneous growth rate 7 
as a function of f is defined through 10''''^ = e^'^ , (or 7(7^) = 27r7(r) log e). Note that the heavy solid 
line for 7 drawn from eq. (72) fits well with various curves computed from the direct numerical 
integrations. The instantaneous growth rates are almost symmetric with respect to their maxima 
near f = 0, as expected. Growth of the modes occurs only when |f| < 18.3, which is in good 
agreement with the results of the direct integrations, demonstrating the validity of the coherent 
wavelet approximation. 

We define a dimensionless amplification magnitude as 



Then, r(f) is an order of magnitude measurement of the amplification of an unstable mode during 
the time interval (—00, f). The total amplification is given by 10'"^°°^. When eq. (75) is substituted, 
the analytic evaluation of the integral in eq. (76) is not an easy task. In view of a shape of 7(f) 
(Fig. 13), we further approximate 7 with a simple form 




(76) 
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where 
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and the termination epoch of growth fc is defined by 

dln.n 



Tc = 



dlnR 




dlni? 



) 4 - «^ - Ql (78b) 



Notice that eq. (77) is valid only if fc is real, that is, only if the condition (73) is satisfied. Prom 
eqs. (76) and (77), the total amplification magnitude is easily found to be 

= rr^' (79) 

which is illustrated with solid contours in Fig. 14. Also shown with dotted contours are the direct 
results from numeral integration of eq. (75), which are in excellent agreement with r(oo). The 
thick contour is the locus of fc = 0, demarcating the stable and unstable regions: the system is 
stable at the right hand side of the thick contour. In the — q^ plane, the total amplification 
tends to be greater as becomes larger and as q^ becomes smaller. This is because the NMRI 
with = acquires maximum instantaneous growth rates at q^ = oo (Fig. 11a) and because the 
shearing time is longer with smaller q^ (cf. eq. [78b]). For comparison, we also include in Fig. 14 
the results from the shearing sheet equations for four parameter sets: (gm, g^) = (0.03, 0.1), (0.1, 
1), (1, 10), and (^2, 100^2), and /?=100 for all cases: these are marked with dots on gm — q^ plane, 
labeled by the respective exact and estimated (in parentheses) amplification magnitudes. Note that 
all of the estimated amplification magnitudes are within 5% of the results of direct shearing sheet 
integrations. This indicates that eq. (79) is an excellent analytic estimate for the amplifications of 
incompressible NMRI modes. 



7.3. Generalized MRI 



Motivated by the success of the coherent wavelet method in finding the solutions of the NMRI 
with purely toroidal background fields, we now generalize both the axisymmetric BH and NMRI in- 
stabilities by considering non-axisymmetric perturbations applied to the rotating medium threaded 
by both vertical and azimuthal magnetic fields. We include the effect of thermal pressure and allow 
the angular velocity O to be a function of i?, but ignore any other radial variations in the initial 
state. We adopt the shearing sheet coordinates as before, and linearize eqs. (1)~(4). After ap- 
plying perturbations in the form of eq. (60), we assume that the perturbations evolve with time 
as e'*'(*)* with the coherent wavelet condition (i.e., d\n^{t) / dt <^ 7(i))- Following the same proce- 
dTirc as §7.2, we obtain the general instantaneous dispersion relation for the MRI (now written in 
dimensional form) 

{2cl + 4)(k • va)'A;2 + {clkl + 4<a(^ + ^z)) + (k • va)A;z^;az^|^ 

"(k-VA)^/c^ + :lid, (80) 
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where k'^ = k^{t) + m? /B? + with the radial wavenumber defined by k^{t) = —mtd^/dR when 

we choose to = —k^{to)/mQ'. When cither m = (axisymmetric case) or v\z = (pure toroidal 
field case), cq. (80) becomes identical respectively with eq. (57) for the BH modes or eq. (72) for 
the NMRI modes. 

From eq. (80) , we obtain the instantaneous instability criteria for the generalized MRI modes 



which are obviously the generalizations of eqs. (59) and (73). It can also be shown that when i = 
90°, both equations (81a) and (81b) become identically v\^{k'^{t) +171"^ /R"^ + k'^) + di}"^ /din R < 0. 
With an O oc rotation profile, eq. (81a) gives the sufficient condition for the instability in an 
extremely cold medium: sin i > {4 — 2a) / (4 — a) . Keplerian flows for example become unstable 
only ii i > 24°, indicating that cold, S^-dominated media are not subject to the generalized non- 
axisymmetric MRI disturbances, just as we found earlier that axisymmetric BH modes are also 
stable in cold flows for small i. We remark that the case with a ^ 2, as potentially possible in 
MHD winds from boundary layers, can just barely satisfy the cold medium instability criterion for 
i ^ 0. Note that unlike the NMRI mode with v^z = 0, maximum growth rates in the Cg 7^ 
case are not achieved as fcz — 00. In fact, high-Zcz or high-m disturbances are efficiently stabilized 
by Alfven and/or MHD waves whenever both poloidal and toroidal fields are present. However 
small they may be, therefore, inclusion of poloidal fields would yield a different result from the case 
with pure toroidal fields (this point was previously noted by Balbus & Hawley 1998). Again the 
stabilizing effect of kinematic shear appears through the time dependence of k^{t), when m 7^ 0. 

Fig. 15a shows how the compressible BH modes are stabilized by azimuthal magnetosonic 
waves. Here, we confine consideration to the radial wavenumber /cr = 0. As Q'm(= v^m/RVt) 
increases, both the growth rates and the ranges in qz{= vp^kz/^) of unstable modes decrease. This 
is because if 7^ 0, azimuthally displaced material feels relatively stronger restoring forces due to 
both thermal and magnetic pressures of the medium as well as stronger tension forces from bent 
field lines. Non-axisymmetric poloidal-ficid BH instability modes become stabilized with increasing 
values of m. When q^^ > y^—dln /din R (= \/3 for a Keplerian rotation), the instability is 
strictly cut off, even when the effect of kinematic shear is not taken into account. 

We remark that the role of temperature of the medium to the BH instability is different 
between axisymmetric (with = and i 7^ 90°; see Fig. 9) and non-axisymmetric (with 7^ 
and i = 90°; see Fig. 15) cases. When i = 90°, as already explained in §7.1, the axisymmetric BH 
instability with = is independent of (3, because only Alfven and sound waves exist and they 
do not interact with each other. When = and i 7^ 90°, magnetic pressure induces vertical 
MHD wave motions which tend to stabilize the system when /3 is small. If /3 S> 1, however, the 
vertical wave motions become nearly acoustic, leaving the toroidal component of perturbed fields 
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unaffected and permitting higher growth rates. When gm 7^ and i = 90°, on the other hand, the 
coupHng of thermal pressure with magnetic pressure occurs through azimuthal MHD motions, and 
the growth rate depends only weakly on p. 

Fig. 15b shows loci of equi-growth rate on the Qz — Qia plane for i = 10° and (3 = 0.01 (dotted 
contours) and [3 = 100 (thin solid contours). For < c/m < 1.17, there exist upper and lower critical 
vertical wavenumbers, Qz^u and q^^\ such that the system is unstable with q^^i < Qz < (?z,u- When 
Qz > Qz,u-, disturbances are stabilized by MHD waves propagating mainly along vertical direction, 
while perturbations with < qz,i approach stable Alfven waves. Each contour has a slope of 
~ — tanz (= —0.18 for i = 10°) at both ends. Note that dotted contours with lower /? are labeled 
with much smaller growth rates than solid ones with higher P, even though they are similar in 
shape. Compared to Fig. 9 or Fig. 15a, this implies that the NMRI instability with a toroidal 
field configuration is more sensitive to temperature than the axisymmetric/non-axisymmetric BH 
instability with poloidal fields. 

When /3 — > DO, from eq. (80) we have the instantaneous growth rates for the generalized MRI 
modes 
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for (k • 'VA^k'^/kz + dO? jdlnR < 0, which is also a generalization of eq. (75). With weak mag- 
netic field strength, one can show from eq. (82) that j~^d'y^/dt ~ —mkji(k-'VA)/{^k^) — as 
(k ■ va)/^ 0. Thus, we see that for the generalized MRIs, the coherent wavelet approach is 
self-consistent in the weak field limit. When the field strength is moderate, on the other hand, we 
obtain d^y'^/dt ~ m^l'kjiCk ■ va)'^ /k'^ . Since 7 ~ ~ (k • va) in this case, the coherent wavelet 
condition is met only when m\Q'\ <C /^r7 (the weak shear limit). Of course, the predominantly 
toroidal-field case that becomes unstable with 3> 1 also satisfies the coherent wavelet condition, 
since k^ becomes arbitrarily large without increasing the (k • VA)-term, as discussed in §7.2.2. 

Comparing eq. (75) with eq. (82), we note that the incompressible MRI can be generalized 
simply by replacing the dimensionless azimuthal wavenumber with (k • va)/^^. Therefore, we 
can write the net amplification magnitude for the generalized incompressible MRI as 
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where the dimensional peak growth rate 70 and the cut-off time of the instability tc are defined by 
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respectively. The net amplification can thus be completely determined by the two parameters: 
the dimensionless wavenumber = (k-VA)/^^ projected in the direction of initial equilibrium 
magnetic fields, and the angle 9 = tan^^ (m / Rkz) of the wavenumber vector with respect to the 
vertical axis. Note that 7o ^ with vanishing qa, while tc ^ oo ss 9 ^ 0, indicating that 
low-m instabilities show higher net amplifications than high-m disturbances as long as qa / 0. 
The total amplification magnitudes, eq. (83), are plotted in Fig. 16 with thin solid contours. For 
comparison, we also plot the numerical results from eqs. (76) and (82) with dotted contours. We 
assume a Keplerian rotation profile. The heavy curve with q\ = 3 cos^ 9 draws the locus of the 
marginal stability. In the limit of a weak magnetic field strength (i.e., 9a ^ 0), it can be shown 
from eqs. (83)~(85) that r(oo) = (2 log e)r2/(/i tan ^), inversely proportional to 9 (for ^ <C 1) 
but independent of qa, as illustrated in Fig. 16. Also shown in Fig. 16 are the results from the 
direct temporal integrations of shearing-sheet equations with /3 = 100 as filled circles (for i = 90°), 
filled triangles (for i = 30°), and open circles (for i = 0°), labeled by the respective exact and 
estimated (in parentheses) amplification magnitudes. These two results agree very well, implying 
that the coherent wavelet approach indeed provides excellent approximations to the solutions for 
amplification of generalized MRIs. 

From Fig. 16, it is apparent that it is the locally near-axisymmetric (in the sense m/Rkz = 
ian9 <C 1) disturbances that experience maximum amplification, with the amplification magnitude 
only weakly dependent on qa = (k • va)/^^ within the unstable regime {qa < 1). The increase 
in amplification factor with Rkz/m predicted from linear theory may in part explain the larger 
amplitudes of power spectra for modes with larger k • z measured from nonlinear simulations of the 
saturated MRI (cf. Hawley, Gammie, &: Balbus 1995). In addition, for the case of pure toroidal 
fields. Fig. 16 suggests only low amplification factors unless is very large, which may help explain 
why Hawley, Gammie, & Balbus (1995) found lower magnetic field saturation amplitudes in cases 
with initial = 0. 



8. Summary and Discussion 

8.1. General Conclusions Based on Linearized Analysis 

Through linear analyses of the ideal MHD equations, we have explored the stability of shearing, 
rotating flows to a wide range of (primarily local) disturbances. The chief motivation for this 
study is to characterize the internal instabilities that could develop in disk winds that emanate 
from an extended region of a differentially rotating protostellar disk around a young star. The 
dynamics of such winds has inspired intensive theoretical effort because they may be responsible for 
observed YSO jets and outflows. In our analysis, we include both results based on generic density, 
magnetic field, and fiow profiles, and results which adopt as initial equilibrium configurations 
the power-law asymptotic solutions of self-confined cylindrically symmetric winds presented by 
Ostriker (1997): p oc R~'^,B^ oc oc and oc v^, <x R~^/'^. For most of our analysis 
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(§§2-6), the flows were assumed to be cold enough that thermal effects can be ignored compared 
with magnetic forces. To make contact with other studies of shear-induced MHD instabilities 
in rotating disks, we also consider stability of specific models which include non-zero thermal 
pressure (§7). For the lowest-order "fundamental" modes, we employ a normal-mode analysis with 
free Lagrangian boundary conditions to find eigenvalues and eigenfunctions of both stable and 
unstable modes (§3). For higher-order modes, we employ three different local techniques to study 
growth of unstable disturbances: In §§4-6, we present numerical and analytic solutions of dispersion 
relations obtained from normal mode analyses in the Rkn ^ 1 WKB limit. These are exact for 
mQ' + ^z^oz ~ disturbances and are valid for a limited time for weakly-shearing circumstances 
where A;r 2> m\i^' kzlRvQ^/vozl (see also §8.2 below). In §7, we employ temporal integrations 
of the shearing-sheet equations to study MRI modes (which are cut off for Rk^ S> 1). We also 
introduce, in §7.2, a "coherent wavelet" formalism which adapts modal analyses for situations where 
shear is considerable (i.e., small Rk^/m); the coherent wavelet analysis is equivalent to a WKB 
approach in the temporal domain. We include a comparison of results from the shearing-sheet and 
coherent-wavelet techniques applied to MRIs, in §7.2.3. 

Applying these techniques we have identified a total of nine different unstable or overstable 
families of disturbances that occur for a wide range of flow parameters: five (FM, BH, ATB, PB, 
and TR) of them are axisymmetric and the other four (NTB, GPB, PR, and NMRI) are non- 
axisymmetric. Table 1 summarizes the properties of these modes. The main general conclusions 
drawn from the analysis in this work can be summarized as follows: 

(1) Systems having a primarily azimuthal magnetic field, for example, disk winds far from 
their source, are susceptible to the fundamental (FM), axisymmetric (ATB) and non-axisymmetric 
(NTB) toroidal buoyancy, non-axisymmetric magnetorotational instability (NMRI) and toroidal 
resonance (TR) modes. Unstable fundamental modes (see §3.2) are concentrated in the central 
parts of jets, and occur in i?(j(,-dominated fiows when the logarithmic gradient of the magnetic field 
is steeper than « —0.75 (cf. eqs. [24] and [25]). Growth rates of unstable FM are comparable to 
inner-wind Alfven frequencies. Long wavelength modes with large amplitudes at large radii are 
all stable, for power-law wind profiles. The TR mode (see §6.1.4) is an overstability, with growth 
suppressed when A;r increases through shear of the vertical velocity, simply becoming oscillatory 
MHD waves. The axisymmetric toroidal buoyancy mode (ATB; §6.1.2) is activated initially by 
the buoyancy force and subsequently by bending poloidal magnetic fields. In geometrical form, it 
is locally similar to the sausage mode of a plasma column confined by toroidal fields, and leads 
to radial mixing. Because growth rates are larger on smaller scales, ATB can contribute to the 
generation of local turbulence in disk winds. The non-axisymmetric toroidal buoyancy mode (NTB; 
see §6.2.1) is much like the Parker instability, but with the centrifugal force replacing the role of 
external gravity. Although the normal-mode analysis for NTB has the largest temporal validity at 
small m/Rk^, the instantaneous growth rate increases with increasing m/Rk^ (cf. eq. [52] and Fig. 
8c). We thus return to the NTB in §8.2, below, applying time-dependent techniques to study the 
Rkji/m <^ 1 limit. Because the NTB is present whenever radial magnetic forces are non-zero, it 
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may be important in promoting radial mixing. Both of the toroidal buoyancy instabilities require 
non-zero magnetic forces in the equilibrium state. The rarefied and cold conditions of disk winds 
do not favor the development of the NMRI (sec §7.2.2 and §7.3). Like the original (poloidal field) 
magnetorotational (BH) instability, NMRI requires dlnQ^/dlnR < 0, but also requires a relatively 
incompressible medium, as is provided by the relatively dense and warm (cg > ^a) conditions in 
an accretion disk. We show the NMRI vanishes in the limit of Cs/va — 0. We further discuss 
perturbations in cold, S^-dominated flows in §§8.2 and 8.3, below. 

(2) Systems having primarily axial magnetic fields, for example, accretion disks or winds very 
near their origin, are susceptible to the Balbus-Hawley (BH), poloidal buoyancy (PB and GPB), 
and poloidal resonance (PR) modes. The well-known axisymmetric Balbus-Hawley instability (BH; 
see §6.1.3 and §7.1) is the most efficient member of the family of magnetorotational instabilities 
(MRIs; see §7.3). It will work to produce channel flows, eventually generating fully-developed MHD 
turbulence through coupling to non-axisymmetric disturbances in the non-linear stage. Driven by 
background magnetic pressure and the centrifugal force, the ax;isymmetric poloidal buoyancy mode 
(PB; see §6.1.1) requires a gradient in the magnetic field strength to be unstable. If the field 
distribution is steep enough, the poloidal buoyancy modes would also work effectively to generate 
radial mixing and turbulence over much smaller scales than the BH instability. Because of their 
overstable characteristics, the impact on the system of poloidal resonance modes (PR; see §6.2.3) 
would be best evaluated with a global rather than local formalism. Configurations with shallower 
background magnetic gradients {q < 1) are also subject to a non-axisymmetric poloidal buoyancy 
instability (GPB; see §6.2.2) which arises in part from geometric effects. 

(3) In distinction to the original, incompressible, axisymmetric BH instability, we found that 
the compressible axisymmetric BH mode is strongly stabilized by the presence of an azimuthal 
magnetic field if the medium has substantially sub-Alfvenic sound speeds. For example, in a cold 
rotating flow with Q oc the axisymmetric BH instability would be completely suppressed 
if the local pitch angle i = tan~^{Bz/ B^) is less than 30° (cf. eq. [47]). In an incompressible 
medium (as provided by a disk with Cg > va), faster sound waves preserve perturbed toroidal 
fields from being dispersed by MHD wave motions, thereby providing a favorable condition for the 
BH instability. When the field configuration is purely poloidal, the compressible BH instability is 
identical with its incompressible counterpart, independent of temperature (cf. eqs. [58] and [59]). 

(4) Even though they share the same instability criterion (cf. eqs. [59b] and [73]), the operating 
mechanisms for the NMRI of purely toroidal B-fields is entirely different from the axisymmetric BH 
instability of primarily poloidal B-fields. In the NMRI (see §7.2), vertical MHD wave motions driven 
by magnetic pressure play an essential role in the feedback loop for induced radial disturbances, 
while the axisymmetric BH instability tends to be stabilized by vertical wave motions. Faster 
sound speeds produce higher growth rates in both instabilities, but for different reasons: in the 
NMRI by activating azimuthal fluid motions preceded by the vertical MHD wave motions; in the 
BH instability by maintaining the perturbed azimuthal fields generated by shear (when Botf, ^ 
0). Because of their non-axisymmetric nature, the NMRI has a transient growth, stabilized by 
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the growth of A;r from kinematic azimuthal shear. For the NMRI mode, we show expUcitly by 

comparison to direct temporal integrations of the shearing sheet equations that the growth rate at 
= can be used to provide a good estimate of the net amphfication magnitude (see §7.2.3). 

(5) The coherent wavelet formalism we develop (§7.2) may be used to compute instability 
criteria and net amplification factors for generalized MRI disturbances with arbitrary magnetic 
field and wavevector orientations (§7.3). Eq. (80) gives the instantaneous dispersion relation for 
generalized MRIs. For strongly compressible flows {cg/vA 0), instability does not occur in S^- 
dominated configurations (cf. eq. [81a]); in this case, flows with an O oc R"^/^ rotation law can 
be unstable only when the magnetic pitch angle i > 24°. Because MHD disk winds generally 
have very small pitch angles, this result has the important implication that such winds will not be 
subject to the development of strong internal turbulence that occurs as a consequence of nonlinear 
MRIs in disks. The absence of MRIs in cold, -B^-dominated winds may be crucial in enabling 
them to propagate over large distances from their sources. High-fez and/or high-m modes of the 
generalized MRI are stabilized by MHD waves, which is a sharp contrast with the NMRI of purely 
toroidal fields in which maximum growth rates are attained at kz — > oo. For incompressible flows, 
the ampliflcation factor for all MRIs can be written analytically in terms of = (k ■ v^)/^ and 
6 = tan~^ (m / Rkz) (eqs. [83]'^[85]); within the unstable regime (q^ < \dln^l/dlnR\^/^, from eq. 
[81b]), the amplification is ~ exp[2ri/(Ktan0)], favoring "locally-axisymmetric" disturbances. 



8.2. Effect of Shear on Dynamical Growth of Buoyancy Instabilities 

Apart from the results of §7 where we adopted the shearing sheet formulation of the dynamical 
equations to study MRIs, the results in this work have been elicited on the basis of the local normal 
mode analyses. As described in §4, these modes may have a limited range of temporal validity, due 
to the effects of background shear. For axisymmetric disturbances with negligible vertical shear (i.e., 
mQ' + kzVQ^ 0), the results presented in §5 and §6 are acceptable for all time; the modes with pure 
imaginary lo will show an exponential growth without interruption over arbitrarily long time until 
nonlinearity sets in. However, for non-axisymmetric disturbances, or for flows with non-negligible 
vertical shear, unstable modes identified in §5 and §6 are not purely growing. As time evolves, the 
differential velocities build up the radial wavenumber through the kinematic shear (cf. eq. [40]), 
which in turn tends to stop the further growth of disturbances. This can be seen directly through 
the suppression of instabilities in the local analysis when k^ is large (cf . Fig. 8) . The characteristic 
time for the wave pattern to change by a fraction e is t = ek^/\mQ' + kzVQ^l; over this interval, 
the disturbance will be amplified by a factor exp(eA:R,Im(a;)/|mr2' -|- fez^^ozD- When k^ 3> m/R, kz, 
Fig. 8 shows that Im(a}) oc k^^, so that the net amplification factor is nearly independent of A;r. 
Since, however, Im(a)) is not larger than ~ jmO' -|- kzv'oJk^^, only order-unity amplification can 
be expected for disturbances which are consistent with the requirements for quasi-steady normal 
mode analysis. 

Because the normal-mode dispersion relations indicate larger values of growth rates when 
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Rkji/m is small (which is however not self-consistent with the WKB treatment), it is desirable to 
extend investigation to allow for Rk^/m small. The coherent wavelet formalism used for the MRI 
in §7 suggests that when m ^ 1 (or Rk^ S> 1 for v'q^ ^ cases), this can be done by regarding kn 
as a time-dependent variable according to eq. (40) and using the asymptotic dispersion relations 
of §6 (i.e., eq. [43] for PB, eq. [45] for ATE, eq. [52] for NTB, and eq. [53] for GPB). To verify 
this argument, we specifically consider the NTB modes (which are one of the chief instabilities 
in 5(^-dominated winds) and compare the results with the shearing sheet temporal integrations. 
For the latter, we set Bz = and integrate eqs. (26)~(31) in time, setting all of the coefficients 
to constant values. The resulting instantaneous growth rates and time evolutions of variables are 
plotted in Fig. 17 as functions of the normalized time r = iO. We omitted the ^R-dependent term 
in eq. (26) in order to remove rapid oscillations arising from a phase mismatch between the density 
and radial velocity; the amplitude evolution is independent of this term. We also neglected the 
vertical velocity shear and and selected q = A;pi,(0) = k^ = 0, R^^ = O.lvj^^, and m = 100. As 
initial conditions, we chose 0.1 for every variable except = 0.01 and integrated the system of the 
linearized equations. Various curves arc computed from direct mmicrical integrations of shearing 
sheet equations, while the heavy solid lines are drawn from the normal mode solution, eqs. (52) 
and (86) (see below), after taking allowance for the time dependence of A;r. The rapid fluctuations 
of the perturbed variables for r < are due to MHD waves with high \k^\, disappearing after 
variables grow substantially. Again, most of growth occurs over a relatively short period of time 
near k^ ~ 0. Note an excellent agreement between the results from two different approaches; we 
have also obtained similar results with integration from other initial conditions. This confirms that 
our normal mode results can also be applied to high-m disturbances if /cr is allowed to vary with 
time. 

Using eq. (52) with k^iit) from eq. (40), we integrate oj over time to estimate the net amplifi- 
cation for the NTB modes 



Xijt) 
Xi(0) 



cxp / Im(c;;) dt 
Jo 



^m + Ji+(V' 



(86) 



where we put v'q^ = kji{0) = and assume a Keplerian rotation. Thus instead of an exponential 
growth, at later time of evolution we have the power-law growth due to the kinematic shear. This 
behavior is distinct from the MRI modes, which are strictly stable for large enough A;r. However, 
the continued local growth of buoyancy perturbations is offset by the role of kinematic shear in 
mixing phases of disturbances. Considering waves on z=constant plane in a square with sides L, the 
maximum averaged contrast in any variable relative to the mean value is {X/ L)xi{t), where A is the 
local wavenumber of the waves. With A ~ kR{t)~^ and using eqs. (6) and (86), the average contrast 
for the NTB modes evolve with time as ~ 

(J7i)V2{l-?)/3-l^ 

vanishing as t — ^ oo for < q < 1. 
Prom eq. (86), amplification factors are essentially scale-free. Although there may be a significant 
growth of the NTB modes on large scales, their dynamical effect on small scales is limited by the 
phase mixing due to shear. 
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In the presence of vertical shear, the evolution of ATB modes is also affected by the kinematic 
growth of kji although they are axisymmetric modes. The net amplification of the ATB modes 
follows a power-law growth as that of the NTB mods does. In fact, eq. (86) with O replaced by 
2vq^/3 gives the temporal behavior of the net amplification for the ATB modes. 



8.3. Discussion of Applications to Protostellar Winds 

In order for a disk wind to overcome the gravitational barrier due to a central object and to 
be centrifugally launched from the surface of a Keplerian disk, the poloidal components of field 
lines should thread the disk at an angle of 30° or more from the axis (Blandford & Payne 1982). 
Once material starts to flow outward along such field lines, it is accelerated primarily in the radial 
direction by the centrifugal force or by the pressure gradient in the toroidal field. Beyond the Alfven 
surface where the local, poloidal component of the flow velocity is equal to that of the Alfven wave 
velocity, the magnetic field is not strong enough to play a role of "a rigid wire", and the inertia 
of gas becomes important, winding up the field lines to be progressively more toroidal. In this 
process, the azimuthal flow velocity decreases below the corotation value. With the increase in the 
azimuthal component of the magnetic field, the associated hoop stress provides the coUimation of 
the outflow and causes the streamlines to bend upward. The radial flow velocity of the outflow is 
still positive, although it decreases gradually, eventually becoming zero at the cylindrical asymptotic 
limit. The power-law solutions (with v^i = -Br = 0) we adopted for many specific cases represent 
the asymptotic limit of each streamline. 

Ostriker (1997) presented self-similar steady solutions for disk winds with cylindrical asymp- 
totics and gave the asymptotic fluid and Alfven speeds and the location of the asymptotic stream- 
lines, characterized by q together with Rx/Ri or Rq/Ri, where Rq, Ra, and i?i denote the radii of 
the footpoint, the Alfven surface, and the asymptote of each given streamline, respectively. Typical 
numerical values for those solutions are = 0.2Qq, v^^/R = 0A2Qq, and v^^/vz = 6 for q = 0.5, 
and Q = O.IQq and vx(f>/R = O.450o for q = 0.9, where Oq is the Keplerian rotation rate at the 
streamline footpoint. 

Using these values we can estimate the growth times of the global fundamental mode and the 
fastest growing (with /cr near 0) toroidal buoyancy modes. The foregoing analysis suggests that 
these disturbances will play the most significant dynamical role, given the ineffectiveness of MRIs 
in cold, S<^-dominated flows. We define the time to grow by V orders of magnitude as tv- For FM, 
ir,FM = r/(|a;| loge), so from Fig. 2 and eq. (25) we have for q = 0.9 

For ATB and NTB, from eqs. (45) and (52), we have for z ^ and q = 0.5 

10^ p /M\-^/V ^0 
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ir,NTB 



3^0 



= 1.7 X 10 



yrs 




10 AU 



i?0 



3/2 



respectively. Here, M is the mass of the central star and Qo.i is the angular speed of a disk at the 
footpoint i?o,i of the innermost streamline of winds. The fact that the growth of the FM by a factor 
10^ occurs within ~ F times the rotation period of the disk at the inner radius, far shorter than 
the lifetime of winds (~ 10^ — 10^ yrs), suggests that the FM mode is dynamically important in 
the evolution of the disk winds. When q is small, the radial turbulent mixing of the wind, caused 
by both axisymmetric and non-axisymmetric toroidal buoyancy modes over a relatively short time, 
is likely to cascade down into arbitrarily smaller scales to dissipate when the microscopic processes 
such as magnetic reconnection are included. The released energy in the dissipation processes may 
heat up the flow, potentially making a significant contribution to the heating of protostellar winds 
and jets. Because the growth rates of buoyancy modes are proportional to the equilibrium magnetic 
forces (cf. eqs. [45] and [52]), winds that have approached a force-free magnetic configuration will 
not be subject to the ATB and NTB instabilities. 

The global fundamental mode affects only the inner region of the disk winds (e.g., the central 
tenth for the model shown in Fig. 3) . The logarithmic density gradient dlnp/d\nR changes relative 
to the equilibrium value by d{pi/ po) / d\nR. Fig. 3 shows that as a consequence of the fundamental 
mode, the very central region becomes more steeply stratified, a surrounding concentric region 
less steeply stratified, and the balance (most of the wind) remains nearly unchanged. Thus, the 
FM tends to enhance jetlike structure in the central parts of winds. In addition, because of their 
tendency to compress interior gas via the FM mechanism, disk winds may help to coUimate any 
interior flows into narrow, fast jets, even when the disk winds themselves have relatively slow motion 
(cf. Ostriker 1997). 

What do the present results imply about the likely radial extent of protostellar winds? First, 
we note that observed optical jets are unlikely to be isolated structures, because if so they would be 
significantly overpressured relative to the ambient medium: Since magneto centrifugal jet models 
typically predict internal Alfven speeds comparable to their flow speeds in the range 150 ~ 400 
km s~^, they have strong internal magnetic pressure Pwind = -^^|/8vr ~ p^v\/2 ~ 2.8 x 10~^ ergs 
cm~^, which is about 6 orders of magnitude greater than the gas pressure of ambient medium, 
-Pext = CgPext ~ 1-3 X lO""*^^ ergs cm~^. Here as reference values we adopted pw = 350 mn cm~^, 
Pext = 200 TTT-H cm~^, = 310 km s~^, and Cg = 0.2 km s~^ with tti-h being the mass of a hydrogen 
atom (Hartigan et al. 1999). Thus the pressure imbalance at the outer boundary of the jet would 
cause either the wind as a whole or only its surface layer to expand until a new balance is attained. 
Based on the results of this paper, if the magnetic field at the base of the wind is stratified less 
steeply than R~^, then perturbations of the outer parts of the wind are stable. As a consequence, 
only the surface layers of such winds would expand in order to achieve a pressure-balanced condition 
with the ambient medium. If, on the other hand, the wind's magnetic field is stratified more steeply 
than at its base, no equilibrium is even possible; the wind would expand as a whole to fill 
the entire 47r steradians, with the inner parts having higher density observable as a narrow optical 
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jet (cf. Shu et al. 1994; Shang et al. 1998). Numerical simulations presently underway (Lee et 
al. 2000) support previous work indicating that protostellar winds with a wide-angle component 
are better able to produce observed molecular outflow structures than purely jetlike winds (see 
also Li & Shu 1996; Ostriker 1997, 1998; Matzner &; McKee 1999) but further studies are required 
to determine just how distributed in angle the wind momentum should be - i.e., to discriminate 
between "fully-expanded" and "surface-expanded" models. Recent observations (see, e.g.. Richer 
et al. 2000) showing a correlation in molecular outflow kinematics with age - with extremely high 
velocity, highly-collimated flows seen only in the youngest sources - may indicate an underlying 
temporal evolution from more-collimated to more-expanded protostellar winds. 

We acknowledge a stimulating report from an anonymous referee, and helpful comments from 
N. Turner and S. Balbus. 

REFERENCES 

Abramowitz, M., & Stegun, L A. 1965, Handbook of Mathematical Functions (New York: Dover) 

Appert, K., Gruber, R., k Vaclavik, J. 1974, Phys. Fluids, 17, 1471 

Appl, S., & Camenzind, M. 1992, A&A, 256, 354 

Bachiller, R. 1996, ARAA, 34, 111 

Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214 

Balbus, S. A., & Hawley, J. F. 1992, ApJ, 400, 610 

Balbus, S. A., & Hawley, J. F. 1998, Rev. Mod. Phys., 70, 1 

Blacs, O. M., & Balbus, S. A. 1994, ApJ, 421, 163 

Blandford, R. D., k Payne, D. G. 1982, MNRAS, 199, 883 

Chandrasekhar, S. 1960, Proc. Natl. Acad. Sci., 46, 253 

Curry, C, & Pudritz, R. E. 1996, MNRAS, 281, 119 

Dubrulle, B., & Knobloch, E. 1993, A&A, 274, 667 

Foglizzo, T., & Tagger, M. 1995, A&A, 301, 293 

Goldreich, P., & Lynden-Bell, D. 1965, MNRAS, 130, 125 

Hardee, P. E., Cooper, M. A., Norman, M. L., & Stone, J. M. 1992, ApJ, 399, 478 
Hartigan, P., Morse, J. A., Tumlinson, J., Raymond, J., & Heathcote, S. 1999, ApJ, 512, 901 



-46- 



Haxtigan, P., Bally, J., Reipurth, B., and Morse, J. A. 2000, in Protostars and Planets IV, ed. V. 
Mannings, A. P. Boss, k, S. S. Russell (Tucson: University of Arizona Press), in press 

Hartmann, L., & MacGregor, K. B. 1982, ApJ, 259, 180 

Hawley, J. F., & Balbus, S. A. 1991, ApJ, 376, 223 

Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1995, ApJ, 440, 742 

Julian, W. H., & Toomre, A. 1966, ApJ, 146, 810 

Konigl, A., k, Pudritz, R. E. 2000, in Protostars and Planets IV, ed. V. Mannings, A. P. Boss, & 
S. S. Russell (Tucson: University of Arizona Press), in press 

Lada, C. J. 1985, ARAA, 23, 267 

Li, Z. Y., & Shu, F. H. 1996, ApJ, 472, 211 

Lee, C. F., Stone, J. M., Ostriker, E. C, k. Mundy, L. G. 2000, in preparation 
Lin, D. C. N., Papaloizou, J. C. B., & Kley, W. 1993, ApJ, 416, 689 
Matzner, C. D., & McKee, C. F. 1999, ApJ, 526, 109 
Miller, K. A., & Stone, J. H. 1999, ApJ, submitted 

Morse, P. M., k Feslibach, H. 1953, Methods of Theoretical Physics (New York: McGraw-Hih) pp 
719-725 

Nagar, N. M., Vogel, S. N., Stone, J. M., k Ostriker, E. C. 1997, ApJ, 482, L195 
Norman, C., k Silk, J. 1980, ApJ, 238, 158 
Ostriker, E. C. 1997, ApJ, 486, 291 

Ostriker, E. C. 1998, in Accretion Processes in Astrophysical Systems, ed. S. Holt, k T. Kallman 
(Woodbury NY:AIP press) p484 

Ouyed, R., k Pudritz, R. E. 1997, ApJ, 484, 794 

Parker, E. N. 1966, ApJ, 145, 811 

Pudritz, R. E., k Norman, C. A. 1986, ApJ, 301, 571 

Rae, I. C., k Roberts, B. 1982, MNRAS, 201, 1171 

Richer, J. S., Shepherd, D. S., Cabrit, S., Bachiller, R., and Churchwell, E. 2000, in Protostars and 
Planets IV, ed. V. Mannings, A. P. Boss, k S. S. Russell (Tucson: University of Arizona 
Press), in press 



-47- 



Reipurth, B., Bally, J., & Devine, D. 1997, AJ, 114, 2708 

Roberts, B. 1985, in Solar System Magnetic Fields, ed. E. R. Priest (Dordrecht: D. Reidel Pub- 
lishing Company) p37 

Rosen, A., Hardee, P. E., Clarke, D. A., Johnson, A. 1999, ApJ, 510, 136 

Ross, D. W., Chen, G. L., & Mahajan, S. M. 1982, Phys. Fluids, 25, 652 

Ryu, D., & Goodman, J. 1992, ApJ, 388, 438 

Shang, E., Shu, F. H., k Glassgold, A. E. 1998, ApJ, 493, L91 

Shu, F. H. 1974, A&A, 33, 55 

Shu, F. H. 1992, The Physics of Astrophysics. II. Gas Dynamics (Mill Valley: Univ. Science Books) 

Shu, F. H., Adams, F. C, & Lizano, S. 1987, ARAA, 25, 23 

Shu, F. H., Lizano, S., Ruden, S., k Najita, J.. 1988, ApJ, 328, L19 

Shu, F. H., Najita, J. R., Ostriker, E. C, Wilkin, F., Ruden, S., & Lizano, S. 1994, ApJ, 429, 781 

Shu, F. H., Najita, J. R., Shang, H., k Li, Z. Y. 2000, in Protostars and Planets IV, ed. V. 
Mannings, A. P. Boss, k S. S. Russell (Tucson: University of Arizona Press), in press 

Stone, J. M., Gammie, C. F., Balbus, S. A., k Hawley, J. F. 2000, in Protostars and Planets IV, 
ed. V. Mannings, A. P. Boss, k S. S. Russell (Tucson: University of Arizona Press), in press 

Terquem, C, k Papaloizou, J. C. B. 1996, MNRAS, 279, 767 

Velikhov, E. 1959, Sov. Phys.-JETP, 36, 995 



This preprint was prepared with the AAS IM^^X macros v5.0. 



-48- 




Fig. 1. — Eigenvalues a = a; -Re /^A (-Re) of stable global fundamental modes for (a) r\ = 0.1 and (b) 
n = 10~^ are plotted as functions of q (i?c is the radius of the outer boundary, Ri is the radius of 
the inner boundary, and r\ = Ri/Re). Each shaded envelope is filled with eigenfrequencies having 
the same n but different i < imax- Since imax is a decreasing function of g, having i^ax = when 
q = 1, the envelopes narrow in width as q increases. Eigenvalues with n = 0.1 are more sensitive to 
q and i than those with n = 10~^. Mode conversion (see text) occurs when a cylindrical wind has 
a narrow thickness (n near 1). Dotted lines in (b) represent asymptotic eigenvalues with n <C 1, 
which are in good agreement with numerical solutions. Also shown are the eigenvalues for the case 
with q = 0.5 and i = as filled circles. The upper boundary of each envelope corresponds to i = 0°. 
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Fig. 2. — Eigenvalues |cj|rj = |ti;|i?i/f A(-Ri) of unstable global fundamental modes. Solid lines 
(n = 10^^) and dashed lines (n = 10~^) are the exact results computed from eqs. (20b) and 
(22). Drawn as dotted lines, the approximate, analytic solutions (eq. [25]) follow the exact values 
fairly well when \a\ is relatively small. The various curves represent different pitch angles of the 
equilibrium magnetic field configuration: i = 0°, 5°, • • •, 35°, 40° from right to left. The uppermost 
thick lines correspond to the maximum pitch angle, imax- Note that |cr|r?^^ is almost independent 
of n. 
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Fig. 3. — Examples of normalized (a) stable and (b) unstable fundamental modes with i = 0° and 
ri = 10~^. We choose q = 0.4, ctq = 2.42, and vir/vo^ = 1 at r = 1 for the stable modes, and 
q = 0.6, |o"|ri^^^ = 0.11, vin/vQ^p = —1 at r = r; for the unstable modes. Stable eigenfunctions 
dominate the outer region, while unstable ones affect only the inner part of the system. 
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Fig. 4. — Normalized frequencies a = CjRo/vx{Ro) of unstable and overstable modes of axisym- 
metric perturbations are plotted against the normalized vertical wavenumber x^, = k^Ro- For all 4 
cases, xji = 10, q = ( = are adopted. In each frame, solid and dotted lines represent imaginary 
and real parts of the solution frequencies, respectively, (a) When i = 0, only the toroidal resonance 
mode (TR) exists, which is overstable with larger real parts, (b) and (c) The TR splits into two 
branches with the inclusion of poloidal fields. Axisymmetric buoyancy modes begin to appear; with 
relatively small i, the axisymmetric toroidal buoyancy mode (ATB) exists, whereas the poloidal 
buoyancy mode (PB) operates with predominant poloidal fields, (d) For a pure poloidal configu- 
ration, only BH and PB modes exist. The BH is dominant at a relatively larger wavelength region 
with a higher growth rate when magnetic field is weak (larger Gr). 
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Fig. 5. — Changes of axisymmetric unstable and overstablc modal growth rates with Gr and i. 
Normalized vertical wavenumber is fixed as Xz = kzRo = 4. A rotation profile O oc is 
assumed and = 10, g = C = are taken. ATB and PB modes need to have smaller Gr to 
be unstable. Note that as Gr increases, corresponding to a weak field regime, only the BH mode 
prevails. Dotted lines in the frames (c) and (d) mark the minimum value of Gr, available for given 
q and i, below which no initial equilibrium exists. For details, see text. 
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Fig. 6. — Frequencies of unstable and overstable modes for non-axisymmetric perturbations with 
171=1 are plotted against the normalized vertical wavenumber. For both frames, = 10, 
g = C = are adopted and an oc rotation profile is assumed. Solid and dotted lines 

represent imaginary and real parts of the normalized wave frequencies, respectively, (a) In a 
toroidal configuration, only TR and NTB exist, with TR being split by the non-axisymmetric 
effect. NTB has a nearly constant growth rate independent of x^- Note that the real part of TR 
increases linearly with Xz, while that of NTB is nearly zero, (b) When i is high, GPB begins to 
appear. PB is almost unchanged by the non-axisymmetry. Like TR, PR is also an overstable mode 
with a larger real part proportional to Xz- 
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Fig. 7. — Effects of i and q on the (a and b) axisymmetric and (c and d) non-axisymmetric buoyancy 
mode frequencies. A rotation profile oc is assumed, and = 10, C = 0, and Gr = 5 are 

adopted. In all four frames, the unstable modes with i < icrit = cos~^ [1 + q) /2 correspond to 
toroidal buoyancy modes, whereas poloidal buoyancy modes have i > icrit- Note that at i = 0, no 
unstable ATB mode exists while its maximum growth rates are achieved at i — > as — > oo. 
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Fig. 8. — Effects of changing = xr/Rq on the unstable/overstable axisymmetric (a and b) and 
non-axisymmetric (c and d) modes. A rotation profile J7 oc is assumed, and Xz = 2 and 

C = for all cases; q = 0.8, Gr = for left frames and g = 0, Gr = 5 for right frames are adopted. 
BH modes cease to exist when ,tr > 3, but all the other modes follow lm{a) ~ at xr ^ 1, 
which is in good agreement with the asymptotic dispersion relations presented in the text. Since 
kinematic shear causes /cr to increase linearly with time, the growth of perturbations would occur 
in a power-law fashion rather than an exponential one at later time of evolution. 
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Fig. 9. — Temperature effect of the compressible BH instability. The abscissa is the normalized 
wavenumber qa = (k-VA)/r2 (= k^VAz/^ for m = 0). A rotation profile O oc R^^^"^ is assumed and 
no radial perturbation is considered (/cr = 0) . Curves with different (3 = cl/v\ show how thermal 
effects modify the BH instability. For i = 90°, the BH growth rate is independent of /?, while (5 has 
a significant impact on the growth rates when i ^ 90°. In a cold MHD limit (/3 = 0), no BH mode 
is expected if z < 30°. 
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Fig. 10. — Schematic diagram showing the development of the NMRI with pure toroidal fields. 
Magnetic fields, thermal pressure, and differential rotation with dVi/dR < all work in cooperation 
to amplify applied disturbances. See text for the detailed explanation. 
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Fig. 11. — Non-axisymmetric magnetorotational instability of toroidal magnetic fields. Here, we 
define = (k-VA)/f^ (= v\^m/RQ = in the text for vaz = 0), qz = VAkz/^ (= VAtpkz/^ for 
1 = 0°), and P = c^/v\^. Rotation with J) oc is assumed. The NMRI instability becomes 

more unstable with (a) higher and (b) higher sound speed. The critical wavenumber q^ is 
independent of /3, but the maximum growth rates are sensitive to temperature for (3 < 10. For 
9m < 1, 7 = gm-\/W(l + f^)- Eventually at ^ = 0, no NMRI is expected. 
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Fig. 12. — An exemplary run of the NMRI shearing sheet equations for q^^ = 0.1, = 1, and 
P = 100. We choose the initial conditions as 6r = 0.01, bz = 0.4, and 0.1 for the other variables, 
and integrate the system of equations from f = —42.44 with which the initial perturbed magnetic 
field is divergence free. Overall evolution can be divided into three stages: initial relaxation phase 
(f < —20), instability phase(|f| < 20), and stable oscillation phase (f > 20). Most of growth occurs 
when |f| is relatively small. Although kinematic growth of the radial wavenumbcr eventually stops 
the further growth of disturbances, the net amplification is tremendous. The coherent wavelet 
solution (eq. [76]) represented by thick solid lines in (a) and (b) is in excellent agreement with the 
shearing sheet results. 
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Fig. 13. — Instantaneous growth rates 7(f) for the NMRI with pure toroidal fields are plotted 
against the orbit f (= ti7/27r) for qin{= vp^^m,/R^) = 0.1, qz{= vp^^k^/Vt) = 1, and (3 = 100. 
The thick curve representing the coherent wavelet solution agrees fairly well with the behaviors 
of various curves for individual variables from the direct numerical integration of the shearing 
sheet equations. The NMRI shows only a temporary growth, but the net amplification is about 9 
orders of magnitude in this example. Kinematic shear increases with time so that eventually it 
suppresses the instability completely at orbit f ~ 18.3. 7(t) is symmetric with respect to f ~ 
and the coherent solution can well be fitted with eq. (77). 
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Fig. 14. — The total amplification magnitude r((X)) of the NMRI with pure toroidal magnetic 
fields is drawn on the qm — Qz plane. Solid contours are computed from the approximate analytic 
estimate, eq. (79), while dotted contours represent the results from direct numerical evaluation of 
eqs. (75) and (76). Four dots correspond to the results of the direct temporal integrations of the 
shearing sheet equations: the adapted parameters are (3 = 100 and ((/m, t/z) = (0.03, 0.1), (0.1, 1), 
(1, 10), and {V2, 100\/2) from lower-left to upper-right. The numbers labeling a dot is the exact 
and estimated (in parentheses) total amplification magnitudes. Note that the analytic estimate 
predicts the true amplification magnitude very closely. r(oo) has a higher value with larger and 
smaller q^- The heavy contour corresponds to the locus of the marginal stability. 
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Fig. 15. — Growth rates of the generalized MRI are drawn as functions of qz = VAk^/^ and 
Qm = VAm/RQ,, (a) for i = 90° and (b) for i = 10°. Wc put /cr = in both frames. In frame (b), 
sohd contours corresponding to (3 = 100 show 7/f^= 0.7,0.6,.. .,0.2, from inside to outside, while 
dotted contours corresponding to /? = 0.01 show j/^= 0.1,0.08,.. .,0.02, from inside to outside. As 
gm increases, both growth rates and unstable ranges of Qz decrease. Eventually, if q^a > 1-73 for 
i = 90° or q^i > 1.17 for i = 10°, the generalized MRI is completely suppressed by MHD wave 
motions. When i = 90°, thermal pressure tends to reduce the growth rates by activating azimuthal 
wave motions if q^ ^ 0; however, the similarity between the f3 = 0, 100 curves in (a) shows its 
effect is not significant. The uppermost thick curve in (b) represents the locus of the marginally 
critical wavenumbers (cf. eq. [81b]) above which no instability can be expected. 
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Fig. 16. — The total amplification magnitude r(oo) of the generalized incompressible MRI is 
drawn as a function of q\ = (h- VA)/f^ and 6 = tan~^ (m/Rkz). Thin solid contours are computed 
from the approximate analytic estimate, cq. (83), while dotted contours represent the results from 
numerical evaluation of eqs. (76) and (82). Both types of contours show r(oo) = 10^, 10^, 10\ 1, 0.1, 
from bottom to top. We adopt a Keplerian rotation profile. The results of the direct temporal 
integrations of the shearing sheet equations with (3 = 100 are shown with different symbols. Open 
circles corresponding to i = 0°, adopted from Fig. 14, are for (q'm, Qz) = VA{m/R,kz)/^ = (0.1, 
1), (1, 10), and (\/2, 100\/2), filled circles corresponding to i = 90° for (gm, gz) = (5 x 10""^, 0.15), 
(0.005, 0.3), and (0.05, 1), and filled triangles with i = 30° for {q„,, g,) = (0.001, 0.2), (0.005, 0.3), 
and (0.1, 1), from left to right. The numbers labeling each symbol is the exact and estimated (in 
parentheses) total amplification magnitudes. r(oo) has a higher value with smaller 9 but nearly 
independent of ^ 1- The heavy curve represents the locus of the marginal stability. 
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Fig. 17. — (a) Instantaneous growth rates and (b) evolutionary behaviors of physical variables 
are displayed for the non-axisymmetric toroidal buoyancy modes with q = 0, Rft = O.lvAtj), and 
m = 100. As initial conditions, we choose 0.1 for all variables except 6r = 0.01. Various curves 
are computed from the direct temporal integrations of the shearing sheet equations. Thick solid 
lines, drawn from the normal mode solution, eq. (52), and its time integration, eq. (86), with 
kR(t) = —mtO,' provide excellent predictions of the numerical results. A Keplerian rotation is 
assumed and vertical shear is neglected. The most significant growth of the NTB modes occurs 
when /cr is small. With increasing /cr, the growth rate is gradually reduced. The rapid oscillations 
for r < is due to MHD waves associated with high |feR|, which are smoothed out as disturbances 
grow. 
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Table 1. Summary of the unstable/overstable mode properties. 



Type 


Geometry of 
Perturbation 


Physical 
Mechanism 


Magnetic field 
Configuration 


Stability 
Character 


FM 


m = kz = 


global mode 


toroidal 


unstable 


ATB 




buoyancy 


toroidal 


unstable 


PB 


m = 


Parker 


poloidal 


unstable 


BH 




MRI 


poloidal 


unstable 


TR 




resonance 


toroidal 


overstable 


NTB 




Parker 


toroidal 


unstable 


PR 


m 7^ 


resonance 


poloidal 


overstable 


GPB 


fez /O 


geometric 


poloidal 


unstable 


NMRI 




MRI 


toroidal 


unstable 



